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1 

Introduction 



Computer control systems are becoming an increasingly competitive factor 
in a wide range of industries. Many products now achieve their competitive 
edge due to the complex functionality provided in algorithms and software. 
As more and more product value is invested in software, there is a strong 
desire to formally verify its correctness. System analysis, which many engi- 
neers may previously have regarded as an academic exercise, is becoming 
instrumental for coping with complexity and guaranteeing correctness of 
advanced software. At the same time, increased performance demands over 
wide operating ranges force control engineers to move from linear to non- 
linear controllers, for which linear analysis techniques often fail. 

Competition also forces faster and more effective product development. 
Today, more and more control designs are based on mathematical process 
models, and their performance is thoroughly tested in simulations before 
full-scale trials. This reduces expensive and time-consuming experimenta- 
tion and tuning on prototype products. Working with mathematical models, 
however, always involves uncertainty: there is always a mismatch between 
what is predicted by mathematical models and what can be observed in re- 
ality. It then becomes important to account for uncertainty in the analysis, in 
order to grant that the results are valid also in reality. 

The amazing advances in computer technology have made high perfor- 
mance computers broadly available. Toaday’s control engineers are skilled 
users of advanced software, and most control designs are today performed 
using sophisticated CAD tools. This makes it very attractive to develop anal- 
ysis and design methods that are based on numerical computations. 

This book presents a computational approach for analysis of nonlinear 
and uncertain systems. We focus on systems with piecewise linear dynam- 
ics and extend some aspects of the celebrated theory for linear systems and 
quadratic criteria to piecewise linear systems and piecewise quadratic crite- 
ria. The analysis tools are easily accessible (as user-friendly computer pro- 
grams) and computationally efficient (relying on convex optimization). 
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Chapter 1. Introduction 

Piecewise Linear Systems 

In this book, we consider piecewise linear systems on the form 

x = AiX + di for x G Xi. 

Here {Xi} is a partition of the state space into operating regimes. The dy- 
namics in each region is described by a system of linear (or, rather, affine) 
differential equations. Piecewise linear systems have a wide applicability 
in a range of engineering sciences. Some of the most common nonlinear 
components encountered in control systems, such as relays and saturations, 
are piecewise linear. Switches also occur naturally, e.g., when a plant op- 
erates in different modes or under physical constraints. Diodes and tran- 
sistors, key components in even the simplest electronic circuits, are natu- 
rally modeled as piecewise linear. Many advanced controllers, notably gain- 
scheduled flight control systems, are based on piecewise linear ideas. The 
construction of a globally valid nonlinear model from locally valid lineariza- 
tions is easy to understand and widely accepted among engineers. 

Some of the first investigations of piecewise linear systems in the control 
literature can be traced back to Andronov who used tools from Poincare 
to investigate oscillations in nonlinear systems [3]. The practical benefits 
of moving from linear to piecewise linear servomechanisms were also no- 
ticed early on [184]. An interesting early attempt to develop a qualitative 
understanding of piecewise linear systems were made by Kalman [111]. He 
considered a saturated system as a series of polyhedral regions in the state 
space, separated by switching boundaries (this is also the view of piecewise 
linear systems that we will adopt in this book). By identifying the singular 
points of the dynamics in the different regions he could then make quali- 
tative statements about the global dynamics. It would take several decades 
before these ideas were refined and developed into more complete analy- 
sis tools. In the meantime, developments on piecewise linear systems ap- 
peared almost exclusively as work on linear systems interconnected with 
static nonlinearities such as relays, saturations and friction. Since these sys- 
tems turned out to be very challenging to analyze, these directions have re- 
mained very active areas of research. Many theoretical results with broad 
applications, such as the work on optimal control [58] absolute stability 
[164], and differential equations with discontinuous right-hand sides [57], 
have their roots in this research. 

It is fair to say that it was the circuit theory community that first rec- 
ognized piecewise linear systems as an interesting system class in its own 
right. Driven by the need for efficient simulation and analysis of large-scale 
circuits with diodes and other piecewise linear elements, a considerable re- 
search effort has focused on efficient representation of piecewise linear sys- 
tems [42, 198, 107]. The research has focused almost exclusively on static 
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problems, while analysis of the more complicated dynamical behaviors has 
remained largely unattended. 

Triggered by the recent increase in the use of switched and hybrid con- 
trollers, two conceptually different approaches to analysis of general piece- 
wise linear dynamical systems have emerged. For discrete-time dynamics, 
Sontag has proposed to exploit properties of affine mappings and polyhe- 
dral sets to formulate decision procedures for answering basic analysis ques- 
tions, see [191]. This approach captures some unique features of discrete- 
time piecewise linear systems but is computationally demanding to imple- 
ment in practice (see the work by Kantner [113] who used similar techniques 
for set-valued simulation of piecewise linear systems). For piecewise linear 
systems with continuous-time dynamics, Pettit has developed a qualitative 
analysis method based on vector field considerations [159]. The approach 
can be seen as a multidimensional extension of phase portrait techniques 
and gives a qualitative picture of the overall dynamics, indicating sliding 
modes, probable limit cycles, and instabilities. More recently, Jirstrand [92] 
and Einarsson [51] have developed alternative analysis methods that com- 
bine local analysis of the systems dynamics within the operating region with 
graph-theoretic ideas to describe the global dynamics. 

This book focuses on quantitative analysis of piecewise linear dynamic 
systems. Stability and gain computations are typical examples. Up until 
very recently, such results existed almost exclusively for particular classes 
of piecewise linear systems (such as linear systems with saturation). The 
research presented in this book has its roots in computational Lyapunov 
theory, which we will describe next. 

Computational Analysis of Dynamical Systems 

Stability analysis of dynamical systems was pioneered by Lyapunov [133, 
134]. The intuition behind his results came from energy considerations: if 
every motion of a system has the property that its energy decreases with 
time, the system must come to rest irrespectively of its initial state. To make 
the argument more rigorous, Lyapunov required that the energy measure 
V (x(t)) of a motion x(t) should be proper in the sense that V (0) = 0, 

V (x) >0 Vx ^ 0 

and V(x) — > oo as ||x|| — > oc. The requirement that V should be decreasing 
along all trajectories of the system x = f(x) takes the form 

V ( x ) = f( x ) < 0 Var^O. 

These are the well-known Lyapunov conditions for asymptotic stability. A 
function V ( x ) that satisfies these conditions is called a Lyapunov function 
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for the system, and can be seen as an abstract measure of the system energy 
Physical insight may sometimes hint at the appropriate energy function, but 
in most cases the choice is far from obvious. To this day, the main obstacle 
in the use of Lyapunov’s method is the nontrivial step of finding an appro- 
priate Lyapunov function. 

In essence, Lyapunov function construction is only simple for linear sys- 
tems. Already Lyapunov showed that a linear system is asymptotically sta- 
ble if and only if it admits a quadratic Lyapunov function V(x) = x T Px. 
The conditions that such a function be proper and decreasing along all mo- 
tions of a continuous-time linear system x = Ax lead to the well-known 
Lyapunov inequalities 

P > 0, A T PpPA < 0. 

With today’s terminology, we would say that these conditions are linear ma- 
trix inequalities (LMIs) in P. Since the inequalities admit an explicit solution, 
this view was not be adopted until almost a century later. Indeed, by pick- 
ing an arbitrary positive definite matrix Q, stability can be assessed from the 
solution P to the system of linear equalities 

A T PpPA = -Q. 

The system is asymptotically stable if and only if P is positive definite. 

Encouraged by these results, several researchers tried to find results of 
similar elegance for nonlinear systems. In particular, much research was fo- 
cused on the absolute stability problem which considers a linear system in- 
terconnected with a static memoryless nonlinearity. The absolute stability 
problem nurtured several important theoretical developments. Two beauti- 
ful examples are the circle criteria [215] and the Popov criteria [163]. These 
results give frequency domain conditions on the transfer function of the lin- 
ear system that are sufficient for existence of certain Lyapunov functions 
for the interconnection. Such frequency domain criteria give valuable in- 
sight and were particularly important before the computer era, since they 
allowed for simple geometrical verification rather than solving difficult ma- 
trix inequalities in the time domain. 

Automatic control went through a drastic change in the 1960’s with the 
advent of state space theory. The development was fueled by demanding ap- 
plications (the space race) , new technology (computers) , and a strong influ- 
ence of mathematics. This led to the development of optimal control, where 
the merit of a control law is often expressed as some integral criteria 

L(x(t), u(t )) dt. 
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Bellman [13] showed that controls u for the system x = f(x,u) that are 
optimal with respect to the above criterion can be characterized in terms 
of solutions V to the Hamilton-Jacobi-Bellman equation 

™ f (^^ /(x)+L(x > w) ) =0 - 

Notice that for the optimal solutions we have V = -L(x, u) which is typi- 
cally negative. Hence V{x) may serve as a Lyapunov function for the closed 
loop system. However, the Hamilton-Jacobi-Bellman equation is notoriously 
hard to solve in general. Many numerical methods have been devised for the 
solution of optimal control problems but they tend to suffer from combina- 
torial explosion. This limits practical applications of optimal control theory 
to systems of low state dimension or to the optimization of trajectories rather 
than feedback laws. 

An important exception is the combination of linear systems and quad- 
ratic criteria. In this case, the dynamics is on the form x = Ax Bu , and the 
criterion takes the form 



x t Qx + 2 u T Cx + u T Ru dt. 

The first solution to this problem was due to Kalman [112] who showed that 
the optimal controller is a linear state feedback. In the early 1970’s, Willems 
gave several equivalent characterizations of the optimal solution [211]. One 
of these characterizations is that there exists a symmetric matrix P = P T 
that satisfies the linear matrix inequality condition 

~A T P + PA + Q PB + C T ~ 

„ >0. 

B t P pc R 

However, the numerical methods at hand made it more attractive to con- 
sider an alternative characterization in terms of an algebraic matrix equa- 
tion which could be solved using numerical linear algebra. Willems also 
showed that many other questions involving quadratic criteria, such as com- 
putations of induced gains, can be characterized by Lyapunov-like functions 
(called storage functions) [212]. For linear systems the existence conditions 
for such functions take the form of linear matrix inequalities. 

A decade later, numerical methods for convex optimization started to get 
widely available. In their 1982 study of the absolute stability problem (now 
extended to multiple nonlinearities) Pyatnitskii and Skorodinskii derived a 
solution in terms of LMIs and gave a numerical algorithm that is guaranteed 
to find a solution when it exists [166]. 




5 




Chapter 1. Introduction 



The early methods for convex optimization had high complexity. A break- 
through came in 1984 when Karmarkar proposed an interior point method 
for linear programming [114]. This method had polynomial complexity and 
worked well in practice. The method was later extended to general con- 
vex programming by Nesterov and Nemirovski [146]. Promoted by efficient 
software [62, 214, 206], strong theoretical developments and excellent tuto- 
rial texts [27, 180] , researchers have started to accept a linear matrix inequal- 
ity condition as a solution of similar value to an analytical result. Moreover, 
linear matrix inequalities have turned out to be convenient for the formula- 
tion of a wide range of important control problems and the interest in these 
methods has soared. 

Since the mid-90’s, there has been a strong activity in trying to extend 
computational results for linear uncertain systems to piecewise linear and 
hybrid systems using non-quadratic Lyapunov functions (see, e.g., [190, 97, 
158, 71, 67]). The combined effort of a large number of researchers has now 
produced a relatively complete set of analysis tools for general piecewise 
linear systems. A key contribution in this book is to show how the search 
for piecewise quadratic Lyapunov functions 





T 


r -i 


X 


p. 


X 


1 




1 



for x e Xi 



for piecewise linear systems can be cast as a convex optimization problem in 
terms of linear matrix inequalities. Based on this result, we extend some as- 
pects of the successful theory of linear systems and quadratic constraints to 
piecewise linear systems and piecewise quadratic constraints. We describe numer- 
ical procedures for assessing stability, computing induced gains and solving 
optimal control problems for piecewise linear systems. These developments 
enable analysis of a large and practically important class of control systems 
that are not easily dealt with using other techniques. The analysis tools are 
easily accessible (as user-friendly computer programs [77]) and computa- 
tionally efficient (relying on convex optimization). 

Philosophy of this Book 

The class of piecewise linear systems studied in this book have nonlinear 
(possibly discontinuous) dynamics and allow switching rules that incorpo- 
rate memory and logic. These systems may exhibit astonishingly complex 
behaviors (including limit cycles, multiple equilibrium points, chaos, etc.) 
and an exact analysis of piecewise linear systems should be expected to be 
hard. In fact, it has been shown that even the problem of establishing sta- 
bility for very simple piecewise linear systems (such as linear systems with 
saturated feedback) is undecidable, see [26]. 
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In light of such results, we cannot hope to find a computationally effi- 
cient (polynomial-time) method for stability analysis of piecewise linear sys- 
tems that is guaranteed to work for all problem instances. Rather, the aim of 
this book is to develop analysis methods that are computationally efficient, 
work for most practical problems, and produce guaranteed results (when 
they work). The methods are based on searching for Lyapunov function in 
a fixed, yet flexible, finite-dimensional Lyapunov function class. The search 
is formulated as a convex optimization problem that can be solved globally 
and efficiently using a wide variety of publically available software. If the 
optimization problem is feasible it returns a Lyapunov function that guar- 
antees a certain system performance. If the optimization problem is infeasi- 
ble, on the other hand, our methods are inconclusive. Infeasibility happens 
when the system does not admit a Lyapunov function within the fixed class 
(this does not necessarily mean that the system is unstable!) or when no such 
Lyapunov function can be found due to conservatism in the computations. 
The approach can be illustrated graphically as in Figure 1.1. 



system 

description 



computer 

program 



— Mg stable 
Mg ‘inconclusive 



Figure 1.1 Computer implementation of approximate stability analysis (from [71]). 

Although the methods are not guaranteed to work for all problem in- 
stances they appear to work very well in practice. They give strict perfor- 
mance improvements over approaches that use quadratic or Lure-type Lya- 
punov functions, and allow analysis of a wide range of systems that cannot 
be analyzed (or are not easily dealt with) using alternative techniques. 

Finally, this book focuses exclusively on continuous-time piecewise lin- 
ear systems. Many results developed in this book have direct analogues for 
systems where the dynamics in each mode is described by affine difference 
equations, see, e.g., [141, 56]. Other aspects of the theory, on the other hand, 
are quite different (see, e.g., [191, 15, 16]). 

About This Manuscript 

This book is based on the author’s Ph.D. thesis [97]. Although this book 
remains faithful to the original at large it has been revised and extended 
in many respects. Several sections have been re-written and many results 
have been improved compared to the original version: issues related to slid- 
ing modes now receive more attention and an alternative analysis method 
based on C 1 Lyapunov functions has been introduced in Chapter 4; Chap- 
ter 7 contains significantly improved analysis procedures for systems with 
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limited regions of attraction as well as a converse Lyapunov theorem for 
smooth nonlinear systems. The chapters on control design and analysis of 
hybrid systems have been extended and new results have been added. This 
book also contains a chapter on convex sets and convex optimization. In 
addition, many chapters include short descriptions of parallel or more re- 
cent work that complement the original content: Chapter 4 shows how Lya- 
punov function computations can be carried out using optimization soft- 
ware that supports equality constraints [72]; Chapter 6 describes a proce- 
dure for quadratic stabilization of piecewise linear systems using globally 
linear control laws [72] and illustrates how design of piecewise linear con- 
trol laws can solved as a bilinear matrix optimization problem (cf. [171]); 
Chapter 8 contains a description of an alternative analysis procedure for hy- 
brid systems that uses Lyapunov functionals [73] and discusses robustness 
issues in hybrid systems [158]. The Comments and References sections of all 
chapters have also been revised to reflect recent progress. 
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Piecewise Linear Modeling 



This book treats analysis and design of piecewise linear control systems. 
In this chapter, we lay the foundation for the analysis by presenting the 
mathematical model on which the subsequent developments will be based. 
We derive an explicit matrix representation of the model and discuss solu- 
tion concepts. We extend modeling techniques for uncertain linear systems 
to piecewise linear systems and show how norm-bound uncertainties and 
smooth nonlinearities can be treated rigorously in this framework. Finally, 
we note that piecewise linear systems enjoy important interconnection prop- 
erties that allow complex systems to be constructed from the interconnection 
of simpler subsystems, and discuss memory-efficient model representations. 



2.1 Model Representation 

A piecewise linear dynamical system is a nonlinear system 



x = f(x,u,t) 
V = 9 &u,t) 



whose right-hand side is a piecewise linear function of its arguments. For 
example, a linear system with saturated input results in system equations 
that are piecewise linear in the input variable u. Linear systems with abrupt 
changes in parameter values (such as jump-linear systems [193] ) are piece- 
wise linear systems in time t. The most common situation, however, is when 
the system equations are piecewise linear in the system state x. Such a model 
can, for example, arise from linearizations of a nonlinear system around dif- 
ferent operating points, or from interconnections of linear systems and static 
piecewise linear components. 

We will understand the term piecewise linear as piecewise linear in the 
system state. With this interpretation, piecewise linear indicates that the 

M. Johansson: Piecewise Linear Control Systems, LNCIS 284, pp. 9-31, 2003. 
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state space can be partitioned into a set of regions, X if such that the dy- 
namics within each region is affine in x 



{ x = AiX + cii + BiU 
y = CiX + c % + DiU 



for xel, 



When written in this way, it is clear that a piecewise linear system has two 
important components: the partition {Xi} of the state space into regions, 
and the equations describing the dynamics within each region. To obtain 
a good understanding of the global dynamics of such systems one needs 
to account for both. Although the term piecewise linear does not impose 
any restriction on the geometry of the regions, such restrictions are often 
necessary in order to arrive at useful results. In this book we restrict our 
attention to polyhedral piecewise linear dynamical systems, where the state space 
is partitioned into convex polyhedra. 



Introductory Examples 

Before stating a precise model definition it is useful to develop a basic un- 
derstanding of piecewise linear systems by studying some examples. We 
will begin with one of the simplest piecewise linear systems:a linear system 
under saturated feedback control. 



Example 2.1— Actuator Saturation in Linear Systems 
Consider a linear system under bounded linear state feedback, 

x = Ax + b sat(v), v = k T x. 

As illustrated in Figure 2.1, the saturation nonlinearity induces a partition 
of the state space into three polyhedral cells corresponding to negative satu- 
ration (Xi), linear operation (X 2 ), and positive saturation (X 3 ) , respectively. 
The dynamics is piecewise linear 



'Ax — b 


x G X\ 




(. A + bk T )x 


x G X2 


(2.1) 


^Ax + b 


x G X3 





In this example, it is natural to let the cells be closed (but unbounded) poly- 
hedra that only share their common boundaries. The presence of offset terms 
makes the dynamics in each cell affine rather than linear in the state x. □ 

Note that the interconnection of a linear system with any other static piece- 
wise linear component (such as a relay, dead-zone, minimum or maximum 
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Figure 2.1 A saturated linear feedback (left) induces a piecewise linear system with 
a polyhedral partition of the state space (right) . 



function, piecewise linear spline, etc.) or linear combinations of such com- 
ponents would also result in a system with piecewise linear dynamics. 

The initial motivation for using piecewise linear models in circuits and 
systems was to approximate nonlinear components in a way that allows 
for tractable computations. This is also a key idea in gain-scheduling ap- 
proaches to modeling and control of dynamic systems. A simple approach 
for constructing a piecewise linear approximation to a smooth function is 
to evaluate the function at a number of grid points and use linear interpo- 
lation between these points to construct the approximant. We illustrate this 
approach by the following example. 

Example 2.2— Approximation of Smooth Systems 

The following differential equations describe a mechanical system with a 

nonlinear spring and damper. 

fxi = /i(x) =* x 2 

\ ( 2 . 2 ) 

= /2W =■ ~X2\X2 \ -Xi(l +X?) 

A piecewise linear approximation of the system equations can be obtained 
by evaluating the right-hand side of (2.2) on the grid shown in Figure 2.2 
(left) and then use linear interpolation between these points. The function 
/ 2 (x) and the piecewise linear approximation / 2 (x) obtained in this way are 
shown in Figure 2.2 (right). □ 

Several other methods for constructing piecewise linear models from data 
have been suggested in the literature, see the Comments and References at 
the end of this chapter. 

With the basic intuition for piecewise linear systems developed in the 
examples, we are now ready to state a more precise mathematical model. 
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Figure 2.2 Partition (left) induced by the grid points marked with x , piecewise linear 
approximation (top right) , and actual nonlinear function (bottom right) . 



Model Definition 

A polyhedral piecewise linear system consists of a subdivision of the state space 
into polyhedral regions and the specification of the dynamics valid within 
each region. We will write the system dynamics as 



(x(t) = Aix(t) + cii + Biu(t) 
= Cix(t) +Ci+ Diu(t) 



for x(t) € Xi iel (2.3) 



Here, x e R n is the continuous state vector, u e R m is the input vector and 
y e R p is the output vector. The matrices A i: B^ c i: Di are constant 

and of compatible dimensions. 

The regions Xi C R n are assumed to be closed (possibly unbounded) n- 
dimensional convex polyhedra which we will call cells . The set of cell indices 
is denoted / and the union of all cells, X = U ieiX it will be referred to as the 
partition . We assume that the cells have disjoint interior so that any two cells 
may only share their common boundary. Depending on the system matrices, 
the right-hand side of (2.3) may or may not be unique on cell boundaries. 

Many results in this book are concerned with the analysis of equilibria. 
Unless stated otherwise, we will assume that the equilibrium point of inter- 
est is located at x = 0. It is then convenient to let /o C / be the set of indices 
for cells that contain the origin and I\ C / be the set of indices for cells that 
do not contain the origin. It is assumed that a* = q = 0 for i e Io. 

Since the cells are closed convex polyhedra, they can be represented as 
the intersection of a finite number of closed halfspaces. In other words, for 
each cell X i} there exists a matrix Gi and a vector gi such that 



— {GiX + gi cz 0 } 



(2.4) 
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2. 1 Model Represen ta tion 



Here, the inequality z y 0 means that all entries of the vector 2 are non- 
negative. The following example demonstrates the notation on the piece- 
wise linear system from Example 2.1. 



Example 2.3— Describing the Saturated Linear System 
Consider the linear system with actuator saturation used in Example 2.1. 
The cells can be represented as in (2.4) with 



G\ — —k T , gi — — 1 , G 2 



k T 

—k T 



9 2 




G3 = k T , g s = 



- 1 , 



The index sets are Iq = {2} and I\ = {1,3}, and from (2.1) we can verify that 
di = Ci = 0 for i € J 0 . □ 



A Matrix Parameterization 

For convenient treatment of affine terms, we define 



x(t) = 



x(t) 

1 



Throughout this book, a bar over a signal vector denotes the augmentation 
of the vector with the unit element 1. Somewhat informally, a bar over a ma- 
trix indicates that it has been modified to be compatible with the augmented 
signal vector, e.g., 



X 




' At 


di 


0 




Olxn 


0 



This allows us to introduce the compact notation 



' A 


Bi ' 


. 


Di 




A 


di 


Bi 


Olxn 


0 


6l xm 


Ci 


Ci 


Di 



(2.5) 

( 2 . 6 ) 



The matrices Si will be called system matrices , and Gi will be called cell iden- 
tifiers . With this notation, the dynamics (2.3) can be re-written as 



x{t) 


— s ■ 


x{t) 


y(t)_ 




u{t)_ 



for {x | GiX y 0} (2.7) 
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and the system (2.3) can be represented by a set of matrix pairs, 
specifying the local dynamics and state space partitioning, respectively 



2.2 Solution Concepts 

A dynamic model can not be fully understood without specifying what we 
mean by a solution to the system equations. One way of defining solutions 
is to specify how to generate the future behavior x(t) of the system from 
any initial state. This approach, closely related to providing a simulation 
algorithm, is easy to pursue for piecewise linear systems with continuous 
right-hand side: in the interior of a cell, integration of the equation (2.3) gives 
unique solutions; when we detect that the solution hits a cell boundary we 
simply determine what cell the state is entering, revise the right-hand side 
of (2.3) accordingly, and continue the integration. 

As we will demonstrate shortly, this approach is not easily extended 
to general piecewise linear systems. For systems with discontinuous right- 
hand side, it it possible that trajectories that enter a cell boundary cannot 
be continued into any of the neighboring regions, or that trajectories can be 
continued into several regions. The following definition of a trajectory al- 
lows us to restrict our attention to the case when the non-smooth dynamics 
does not create any problem for our analysis. 

Definition 2.1— Trajectory 

Let x(t) G U ie jXi be an absolutely continuous function. We say that x{t) 
is a trajectory of the system (2.3) on [to,tf\ if, for almost all t G [to,tf], the 
equation x(t) = Aix(t) + a* + Biu(t) holds for all i with x(t) G X it □ 

Note that this definition of a trajectory does not require that the system has 
continuous right-hand side. If x(t) at some time tk passes through a cell 
boundary where the vector fields in the neighboring regions do not match, 
the model (2.3) does not provide a well-defined time derivative. If there are 
sufficiently few crossing times (technically speaking, if the set of times of 
discontinuity have measure zero [175]), however, these time instants can be 
removed without disqualifying x(t) from being a trajectory, see Figure 2.3. 
Trajectories are allowed to remain on cell boundaries only if the vector field 
in the neighboring regions match. 

As illustrated by the following example, trajectories in the sense of Defi- 
nition 2.1 might be non-unique or not even exist. 
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linear system. The times , marked with dashed lines in the time plot, are the times 
where x(t) G Ii n I 2 and the time derivative of (2.3) is not defined. As the number 
of switches is countable, x(t) still qualifies as a trajectory. 



Example 2.4— Non-existence & Non-uniqueness of Trajectories 
Consider the piecewise linear system 1 



x\ = — 2x\ — 2x 2 sgn(xi) 
x 2 = x 2 + 4xiSgn(xi) 



( 2 . 8 ) 



The state-space partition and the vector field in the interior of the two cells is 
shown in Figure 2.4. We see that trajectories with initial values x(0) e = 



? / / / , A 

/ / / / ✓ k 
/ / / / /, 
f / / // , 
/*///, . 


k vn\W 

cv \ \ \ \ 
. N \ \ N 
. 


* / t i ; r 

f t t i . . 
t t t . , . 

1 ; M 


~ l cvA ' 

. > i \ \ ^ 

. , t t \ 1 



Figure 2.4 Simple piecewise linear system illustrates that trajectories in the sense of 
Definition 2.1 need not be unique or even exist. 

{x | x\ = 0 A x 2 < 0} are not unique since they can be continued into either 
cell. Trajectories that reach the manifold S ^ = {x | x\ = 0 A x 2 > 0}, on the 
other hand, cannot be continued into any of the regions (the vector fields in 
both cells point toward S f). In fact, the only trajectory for this system on 
[0, oo ) is the trivial x = 0. □ 

^ere, sgn(-) denotes the set-valued function sgn(;Q = 1 if z > 0 and sgn(z) = — 1 if z < 0. 
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Non-uniquness and non-existence of trajectories are discomforting from a 
model validation point-of-view. These phenomena indicate that important 
characteristics of the underlying physical system have been neglected and 
raise questions about how to interpret the model equations (see the Com- 
ments and Reference section at the end of this chapter for more details). 
From an analysis point-of-view, however, the main obstacle will be the cases 
when no continuation of a trajectory in the sense of Definition 2. 1 is possible. 
The following definition allow us to single out such situations. 

Definition 2.2— Attractive Sliding Mode 

The system (2.3) is said to have an attractive sliding mode at x s if there exists 
a trajectory with final state x s but no trajectory with initial state x s . □ 

For sake of clarity, we will present the main results in this book for systems 
without attractive sliding modes. Solution concepts for systems with slid- 
ing modes are discussed below and a method for detecting sliding modes 
will be described in Section 3.3. The general analysis procedures will be de- 
veloped for systems without attractive sliding modes and a technique for 
extending these results to systems with attractive sliding modes is given in 
Section 4.8. 

Generalized Solutions 

It is still possible to define meaningful solution concepts for systems with 
attractive sliding modes if we consider functions x(t) that remain on the cell 
boundaries for some time interval. The most well-known solution concepts 
are due to Filippov [57] and Utkin [197]. Typically, solutions are defined by 
some limiting process, such as introducing a small hysteresis around cell 
boundaries and letting the hysteresis parameter tend to zero; see Figure 2.5. 
In this way, the behavior on the surface of discontinuity is defined by av- 
eraging the dynamics in the neighboring regions. The following definition 
specializes the approach of Filippov to piecewise linear systems. 

Definition 2.3— Filippov Solution 

Let x(t) g U ieiXi be an absolutely continuous function. We say that x(t) is 
a Filippov solution of the system (2.3) on [to,t /] if, for almost all t e [to,tf], 

x(t) € co {A k x{t) + a k + B k u(t)} 
keK(t) 

where K(t) is the set of cell indices k such that x(t) G X k . □ 

In the above definition, co denotes convex closure. In other words, x(t) is a 
Filippov solution of (2.3) if and only if there exist positive scalars a k (t) > 0 



16 




22 Solution Concepts 



f f / * ^ n \\N 
'? / / * ^w\N 

? / / S (j\ S\\N 

/ / / ^/V N \ N 
f / t f . A \ \ \ 
? t tf * .. \ \ \ 

* f A - ■ ■ ' r ' 

* t >. i \ * 



/ / / s J 




. v S \\ 


// / r . 




,S\W 


// / ✓ . 




w v\\S 


/ / / / . 




, s \ \ \ 


f / / , s 




. X \ \ S 


ft!/. 




♦ > v \ N 


f f . 




. 1 v \ S 


f r 4 . < 







i r. , »i.r . t l 



* 


v s\W 


/ / / , s 


v^N\\ 


/ / / s , 


X s N \ V 


f / / * J 


L \ \ \ \ 


? f r ,/ 


. N \ \ N 


t f t / . 


* ^ s \ N 


f t /, . 


. , \ \ > 


f t ^ .6 - 
' / f ■ J 


H) , v \ s 

. i 1 ^ 



Figure 2.5 Filippov solutions may remain on cell boundaries following dynamics 
generated by averaging the dynamics in the neighboring regions. 



with a k(t) = 1 such that 

Mt) = m ak O L 4 * 3 ^) + a k + B k u(t)} (2.9) 

keK(t) 

for almost all t G [to, t/]. In the interior of cells, the differential inclusion con- 
tains only one element and Filippov solutions coincide with the trajectory 
concept defined above. Time functions x{t) that remain on a cell boundary 
for some time are also accpeted as Filippov solutions if they can be con- 
structed from some convex combination of the dynamics in the neighboring 
regions. As illustrated in the following example, the additional information 
that x(t) remains on the boundary can be used to narrow down the admis- 
sible behaviors and sometimes even derive unique dynamics on the sliding 
mode. 

Example 2.5— Sliding Mode Dynamics on Simple Boundaries 
Consider again the piecewise linear system (2.8) on the set S± = {x | x\ = 
0 A X 2 > 0}. Filippov’s definition gives 

x(t) G co {Aix(t), A2X2 (£)} = ai(t)Aix(t) + a 2 (t)A 2 x(t) 

where the weights satisfy a\ (t) > 0, <22 (t) > 0 and ai(t) + <*2 (t) = 1. How- 
ever, not all such weights give rise to solutions that stay on In order for 
x(t) to stay on S we must have x\ (t) = 0, i.e., 

an (t) [-2xi (t) - 2 x 2 (t)\ + a 2 (t) [—2xi (7) + 2x 2 (7)] = 0 x(t) G S 

Combining this requirement with (t) + a 2 (7) = E we find 

Oil (t) = OL 2 (t) = 1/2 
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Hence, the unique dynamics on the sliding mode is 

xi(t)=Q x 2 (t) = x 2 (t) x(t)eSi 



□ 

On simple boundaries of piecewise linear systems, it is always possible to 
derive a unique (but not necessarily linear) sliding mode dynamics by bal- 
ancing the vector fields as above. However, sliding modes do not only occur 
on simple cell boundaries but can also arise on the intersection of multiple 
boundaries. As demonstrated by the following example, the sliding mode 
dynamics on intersecting boundaries is in general not unique. 

Example 2.6— Non-Unique Sliding Mode Dynamics 
Consider the piecewise linear system 



xi =x 2 -sgn(xi) 

X 2 =X 3 -sgn(x 2 ) 

x 3 = -2xi - 4x 2 - 4x 3 - x 3 sgn(x 2 )sgn(xi + 1) 



The two first state equations reveal that sliding occurs on the sets 

Si = {x | Xi = 0 A |x 2 | < 1} S 2 = {x | x 2 = 0 A |x 3 | < 1} 

As indicated in the simulation shown in Figure 2.6, sliding can also occur on 
the intersection of these surfaces, 

<Si 2 = {x | xi = 0 A x 2 = 0 A |x 3 | < 1} 



To study the sliding motion on <S 12 in more detail, note that the sign 
functions induce a partition of the state space into four cells, being the four 
quadrants in the xi - x 2 plane (see Figure 2.6, right). The vector fields in 
the four regions point toward the set <S 12 so the sliding mode is attractive 
(solutions reaching <S 12 cannot be continued into any cell). Filippov solutions 
satisfy 



4 

x = ^2 a k {AkX + a k } (2.10) 

k= 1 

Here, the weights a k should be non-negative, sum to one 



Ul + OL 2 + CK 3 + CK4 =1 
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Figure 2.6 A Filippov solution enters and slides into the set S 12 where a new 
sliding motion can be sustained (left) . The rightmost figure shows the vector fields on 
the subspace {# | #3 = 0}: the vector fields in the four cells point toward S12 . 



and guarantee sliding on S 12 , i.e., that x\ = x 2 = 0 when aq = x 2 = 0: 

X \ = 0 => — Gq + U 2 + CK 3 — U 4 =0 

x 2 = 0 => x 3 - ai - a 2 + Ois + a 4 =0 

Since there are four vector fields to balance but only three constraints, one 
might suspect that the dynamics in the sliding mode is not unique. Indeed, 
it is straightforward to verify that one solution is given by 

1 + sgn(x 3 ) l+x 3 x 3 1 

ai = , a 2 = — cti, u 3 = — — + ai, a 4 = - - oq 

while another solution is 

aq = a.2 ol<i = aq aq = a 4 ol 4 = a 3 



In fact, these weights define the Filippov solutions with the minimal and 
maximal velocity for x 3 on <Si 2 . The time responses shown in Figure 2.7 re- 
veal that the two definitions give rise to drastically different sliding mode 
dynamics, with settling times differing by at least a factor of four. □ 

Non-uniqueness of the sliding mode dynamics indicates that some impor- 
tant characteristics of the underlying physical system has been neglected, 
and that one needs to review the assumptions made when deriving the 
model. A unique sliding dynamics can sometimes be recovered by model- 
ing salient features of the switching mechanisms, such as the relative speed 
of relay components (see, e.g., [136]). 
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Figure 2.7 Filippov’s definition yields non-unique velocities along the sliding mode. 



The fact that sliding modes can occur on intersections of multiple bound- 
aries makes it hard to verify that a given model is free from attractive slid- 
ing modes. Not only is there an exponential number of surfaces that need 
to be checked for potential sliding modes, but it is also much harder to 
verify that sliding modes on intersecting surfaces are attractive (i.e., that 
solutions starting close to an intersection converge to it). In essence, one 
then needs to do a complete stability analysis of the potential sliding mani- 
fold, which quickly becomes prohibitively expensive. The following exam- 
ple, taken from [201, 197], illustrates this point. 



Example 2.7— Establishing Attractivity of Sliding Modes 
Consider the following piecewise linear system 



xi = — sgn(xi) + 2sgn(x 2 ) 
x 2 = -2sgn(xi) — sgn(x 2 ) 



The sign-functions induce a partition into four cells, each being one of the 
quadrants in R 2 . The vector field is constant in each cell and it is strongly 
suggested that solutions should behave as shown in Figure 2.8: solutions 
spiral toward the origin, which is reached in finite time, and remain there. 
Since the vector fields within the four cells do not match, the origin is a 
sliding mode, and the Fillipov solution x(t) = 0 can be constructed by giving 
the system matrices in the four cells equal weights in (2.9). 

It is far from trivial to prove that the origin is attractive. Vector field con- 
siderations are not easily used since no vector field points towards the ori- 
gin. For this particular system, convergence to the origin can be proved by 
noting that d/dt(\x\(t)\ + \x 2 (t)\) = -2, which means that solutions with 
initial state (xi( 0 ),X 2 ( 0 )) cannot stay away from the origin for longer than 
(|xi(0)| + |a ?2 (0)|)/2 units of time. □ 
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Figure 2.8 Piecewise linear system with non-traversal sliding. 



2.3 Uncertainty Models 

Uncertainty and robustness are central themes in modeling and analysis of 
feedback systems. One of the most important reasons for using feedback 
is to guarantee that system specifications are met despite variations in sys- 
tem components and exogenous disturbances. Furthermore, since there is 
always a mismatch between the model used in the mathematical analysis 
and the actual physical system, it is important to account for this uncer- 
tainty to ensure that the results derived from the model will hold in reality. 
In this section, we will extend the standard uncertainty models for linear 
uncertain systems (see, e.g, [220]) to systems that are piecewise linear. These 
uncertainty models will enable us to extend analysis procedures for piece- 
wise linear systems to yield rigorous results for smooth nonlinear systems. 

To verify robustness we have to specify the sets of admissible uncertain- 
ties and disturbances. We will consider two main classes of uncertainties. 
The first deals with systems on the form 



X = f{x) 



and considers the effects of uncertainties in the function f(x). This situation 
may occur when f(x) is a piecewise linear approximation of some smooth 
function. If the uncertainty is due to unknown or time- varying parameters, 
this is usually called parametric uncertainty. The second class of uncertainty 
descriptions considers systems on the form 

x = f(x,y) 

y = g(x,y) 

where g{x , y) is uncertain or lacks a description with appropriate structure. 
This type of uncertainty is usually called dynamic uncertainty , and may occur 
when y represents an exogenous disturbance or a neglected component. 
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Piecewise Linear Differential Inclusions 

One way to embed more general nonlinear systems in the piecewise linear 
framework is to allow systems with time-varying system matrices 



x(t) = Ai(t)x(t) + di{t) + Bi(t)u(t) 
y(t) = Ci(t)x(t) + Ci{t) + Di(t)u(t) 



for x G Xi 



We will consider the case when the system matrices Si for each cell can be 
written as a convex combination of matrices S'/, ... , S/C In other words, we 
assume that for every t there exist scalars a k (t) > 0 with J2 k a k(t) = 1 such 
that Si(t) can be written as 



K 

S i (t) = J2*k(t)S*. (2.11) 

k = 1 

We will then consider the family of models obtained by considering all ad- 
missible a k (t). For notational convenience, we will to each cell Xi associate 
an index set K(i) that specifies the matrices that are used in the inclusion. 
We will then write (2. 1 1) as 



(2A2) 

We will call these models piecewise linear differential inclusions , pwLDIs. 

An absolutely continuous function x(t) is called a solution of the inclu- 
sion on [to,tf] if, for almost all t e [to,tf], it satisfies 



x(t) 

y(t) 



G CO 
k£K(i) 





for x(t) eXi (2.13) 



Linear differential inclusions have been used to model parametric uncer- 
tainty in linear systems. One of the most well-known examples is the sector 
conditions that have been used in the work on absolute stability [132, 215, 
163]. In this context, the extension to piecewise linear differential inclusions 
allows us to use piecewise linear sector bounds to embed smooth nonlinearities 
into the piecewise linear framework. 



Example 2.8— Sector Bound Nonlinearity 

Consider an integrator in a negative feedback loop with static nonlinearity 
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x(t) = -<p(x(t)) 
y(t) = x(t) 




2. 3 Uncertainty Models 





Figure 2.9 Analysis of smooth nonlinear system (left) via piecewise linear sector 
bounds on the nonlinearity (right) . 



Assume that the nonlinearity <p(x(t)) can be bounded by “piecewise linear 
sectors”, see Figure 2.9. In other words, assume that there exist vectors k 
and Hi, describing upper and lower bounds respectively, such that 

<p(x(t)) G co {T[x(t), ujx(t)} for x G Xi. 

The closed loop system can then be described by the piecewise linear differ- 
ential inclusion 



s^emis^sr} 



with 
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Notice that the integrator system could be replaced with a general piecewise 
linear system, and all results would follow similarly. □ 

Methods for obtaining upper and lower piecewise linear bounds on scalar 
functions can be found in [129] and a procedure for identifying tight up- 
per and lower bounds from data has been proposed in [107]. We will use 
the piecewise linear sector bounds to analyse smooth nonlinear systems in 
Sections 4.7 and fuzzy control systems in Section 7.4. 

Norm-Bound Approximation Errors 

One problem with uncertainty descriptions in terms of differential inclu- 
sions is that the analysis conditions must consider every extreme dynamics 
in each region. A careless application of pwLDIs in the modeling phase can 
generate a large number of extreme systems which may render the analysis 
computations intractable. 
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If the piecewise linear system is obtained by approximating a smooth 
system, it is natural to use a norm bound on the error 

||/(x) - AiX - ai\\ < ti | x | for x € X i7 i e I. 

between the right-hand side of x = f(x) and its piecewise linear approxima- 
tion. Now, rather than specifying several extreme dynamics for each cell we 
only have to provide a single norm bound e*. Moreover, if f(x) is smooth, the 
bounds on the approximation error will decrease as the partition is refined. 
In Chapter 7, we will use this observation to establish a converse theorem 
for analysis of smooth nonlinear systems using piecewise linear techniques. 

Dynamic Uncertainties and Dissipation Inequalities 

The standard approach to account for dynamic uncertainties is to consider 
norm-bound uncertainties in a feedback interconnection, as shown in Fig- 
ure 2.10. In robust control literature, the nominal system E is assumed to be 





A 












W 


2 
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Figure 2.10 Piecewise linear system with uncertainty feedback. 



linear time invariant, while system nonlinearities and uncertainties are con- 
fined to the uncertainty block A. In contrast, we will allow E to be piecewise 
linear. This allows us to choose whether to describe system nonlinearities 
explicitly in the piecewise linear subsystem or to model them as uncertain- 
ties in the A subsystem. This additional freedom can be used to trade off 
computational complexity against conservatism in the analysis. 

The operator A that specifies the feedback w = Az may be linear time 
varying or nonlinear, but is assumed to satisfy the dissipation inequality 




z ( s ) 
w(s) 



M 



z(s) 

w(s) 



ds > 0 



for all t > 0 



(2.14) 



for some real symmetric matrix M . This includes, for example, passive com- 
ponents and elements with bounded £2 -induced gain. After establishing 
gain and passivity properties of E and A using, for example, the techniques 
in Chapter 5, we may try to use small gain or passivity results to establish 
stability of the closed loop system [116]. 
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A Modeling Trade-Off: Uncertainty versus Complexity 

The piecewise linear sector bounds allow stepwise refinements of a global 
sector bounds to improve a nonlinearity description. In the end, such refine- 
ments allow arbitrarily tight inclusions of any continuous function. How- 
ever, each such refinement comes at the price of increased memory require- 
ments for the model representation and increased computational cost for 
the analysis. It is thus natural to look for the simplest model that gives a suf- 
ficiently accurate answer in the analysis. The following example illustrates 
this idea. 

Example 2.9— Piecewise Linear Modeling and Complexity 
In Chapter 5 we will encounter a system on the form 

x(t) = Ax(t) — B(f(Cx(t)) 



where the nonlinearity <g(-) is a spring with the piecewise linear character- 
istic shown in Figure 2.11 (left). Nonlinear spring arrangements of this type 
can be found in engine control systems (see, e.g., [120]). The exact descrip- 
tion has many segments and results in relatively demanding analysis com- 
putations. It is therefore attractive to base the analysis on an approximate 
model that requires less computations. The piecewise linear sector bounds 
give many possibilities. For the specific example in Chapter 5, the piecewise 
linear sector bounds shown in Figure 2.11 (middle) can be used to asses sta- 
bility while an analysis based on global sector bounds (right) fails. □ 




Figure 2.11 From exact description to global sector bounds. Piecewise linear sector 
bounds allow us to trade precision of the model against simplicity of the computations. 



Although the principle of refinements using piecewise linear sector bounds 
has been illustrated on a scalar nonlinearity the methods applies also to mul- 
tivariable nonlinearities. 
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2.4 Modularity and Interconnections 



Modularity and structure-preserving interconnections are attractive features 
in modeling and analysis of dynamic systems. Modularity allows complex 
systems to be represented as the interconnection of simpler subsystems. Im- 
portant model components can be stored in a library and recalled and in- 
terconnected as needed. Structure-preserving interconnections are also very 
useful since they grant that the interconnected system shares important struc- 
tural properties with its components. For example, series, parallel and feed- 
back interconnections of linear systems are themselves linear systems. This 
allows the full systems and its components to be analyzed using the same 
tools. We can choose whether to analyze the full system directly, or to first 
analyze its subcomponents and then invoke interconnection results, such as 
small gain and passivity theorems. 




Figure 2.12 Series, parallel, and feedback interconnection of dynamical systems. 

The following proposition states that the most common interconnection 
structures preserve the piecewise linearity of its components (cf. [198, 159]). 

Proposition 2.1— Structure-Preserving Interconnections 
Series, parallel and feedback interconnections (without algebraic loops) of 
polyhedral piecewise linear systems are themselves polyhedral piecewise 
linear systems. Moreover, a matrix representation {S^Gi} ie i for the total 
system is obtained directly from the matrix representations of its compo- 
nents. □ 

Proof: The result follows by direct computations, see Section B.l. 

The proof of Proposition 2.1 reveals that the interconnection of two piece- 
wise linear systems results in a combinatorial growth in the number of cells. 
Consider the interconnection of the system £i with state vector x and par- 
tition {Xi} ie i and the system £2 with state vector 2 and partition {Zj}j e j. 
The partition of the interconnected system is obtained by considering all 
combinations of (i, j) such that xel* and z £ Zj. Hence, if £1 has a parti- 
tion of pi cells and £2 has a partition ofp 2 cells, the interconnected system 
may have pi x P 2 cells. This illustrates the usefulness of an input-output 
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Figure 2.13 Series connection of two piecewise linear systems is piecewise linear 
(left) . The interconnected system has the partition shown to the right. 



analysis for piecewise linear systems: if analysis of the interconnected sys- 
tem is too computationally demanding, we can try to analyze the simpler 
subsystems and apply small gain or passivity arguments. 

Proposition 2.1 allows descriptions of standard piecewise linear compo- 
nents, such as relays and saturations, to be stored in a library and recalled 
when needed. In this context, it is useful to allow partitioning in the product 
space of the input space and the state space, i.e., to let 
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Although this route will not be pursued here we note that a similar result to 
Proposition 2.1 holds also in this case. The following example illustrates the 
interconnection properties of piecewise linear systems. 

Example 2.10— Interconnection Preserves Structure 
Consider the series connection of two piecewise linear systems 

x = — sat(x) + u 
V = -sgn(y) - x 

illustrated in Figure 2.13 (left). The interconnection is itself piecewise lin- 
ear in accordance with Proposition 2.1. The individual subsystems have two 
and three cells, respectively, while the interconnected system has six. Al- 
though the combinatorial growth is not so pronounced in this example, it 
can be significant in more complex systems. □ 

Unfortunately, the modular approach does not extend directly to the case 
when the systems are uncertain in the pwLDI-sense, since this may intro- 
duce uncertainty in the state space partitioning. Analysis of piecewise linear 
systems with uncertain partitions will be discussed in Chapter 8. 



27 







Chapter 2. Piecewise Linear Modeling 

2.5 Piecewise Linear Function Representations 

We have proposed to represent piecewise linear systems using the matrix 
parameterization (2.7). This representation is easy to understand and con- 
venient to use for analysis. However, if the piecewise linear system (2.3) has 
continuous right-hand side, the representation contains a lot of redundant 
information and is not very memory efficient. 

Modeling and simulation of piecewise linear systems has attracted a 
large interest in the circuits and systems community during the last decades. 
Driven by the need to simulate large-scale circuits with piecewise linear 
components, a large research effort has been devoted to deriving memory 
efficient representations for piecewise linear systems [42, 126, 107]. More 
compact descriptions than the matrix representation (2.7) can be obtained 
when the vector field of the system is continuous across cell boundaries. To 
see this, consider the situation in Figure 2.14. To obtain continuity on the cell 




Figure 2.14 Continuity of vector fields allows parameter savings. 



boundary dXij 



{x | gfjX + hij = 0}, the system matrices must satisfy 

Aj = Ai + Cij gjj 



for some e R n . Since the boundary equations of the cells re-appear in the 
description of feasible changes in the mapping, there is a certain redundancy 
in the data given by the simple model (2.7). The argument above indicates 
that it should be sufficient to store one linear system description, the bound- 
ary equations, and the update vectors c^. 

The first compact parameterization of that appeared was the canonical 
piecewise linear function description introduced in [42]. For piecewise linear 
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dynamic systems, it takes the form 

p 

x = Ax-\-a -\-^2ci\gfx-\- hi\. (2.15) 

i= 1 

This representation stores only a single affine system description, the bound- 
ary hyperplane equations and the update vectors. It has also eliminated the 
need for explicit storage of cell identifiers and no cell identification is nec- 
essary to evaluate the mapping. The representation (2.15) is very efficient 
compared to the simple matrix parameterization (2.7), but it can only rep- 
resent a subset of the continuous piecewise linear mappings (cf. [109, 115]). 
To overcome this problems, various higher order basis function expressions 
have been suggested. These are much more complicated than the simple 
model (2.15), but can be put in the general form 

q 

x = Ax + a + Ci(pi(x) 

i= 1 

Here, the cpi(x) are piecewise linear functions constructed from nested abso- 
lute or maximum functions, see [109, 70, 115, 107]. 

An alternative formulation, which is closely related to the matrix repre- 
sentation (2.7), is the implicit piecewise linear function description [198, 115] 

f x = Ax + a + Bu 

< i = Gx + g + Cu (2.16) 

^0 = u T i u >z 0, i h 0 

This representation was derived from a static linear network where some 
ports have been terminated by negative ideal diodes. The variables i and 
u correspond to currents and voltages respectively, and the last equation 
of (2.16) describes the characteristic of an ideal negative diode. Vectors u 
and i that satisfy this equation are called complementary. To see the close 
connection to our matrix parameterization, consider the region where u = 0. 
Then, the above model reduces to 



x = Ax + a for x such that Gx + g y 0. 

Given x, an evaluation of the mapping needs the associated diode volt- 
ages u. These can be computed by solving the linear complementary prob- 
lem of finding i and u that satisfy the last two equations of (2.16). In princi- 
ple, a solution to this problem can be obtained through a sequence of piv- 
oting operations around the C-matrix. Each such pivoting step forces one 
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entry of the vectors u or i to zero while the corresponding entry in the other 
vector is allowed to be non-zero. The second equation in (2.16) implies that 
the zero entries of i and u force the corresponding entries of the other vec- 
tor to be affine expressions of x. Hence, every configuration of u and i into 
zero and non-zero entries define a polyhedron in the state space (via the 
conditions i >z 0 and u y 0). In each polyhedron, the affine expressions for 
u describe the changes in the local dynamics via the first equation of (2.16). 
In this way, the matrix C encodes the changes in the cell descriptions while 
the matrix B encodes the changes in the affine dynamics. For more details, 
see [198,38, 159, 107]. 

One drawback with the implicit model (2.16) is that a solution to the lin- 
ear complementary problem may require a number of pivoting operations 
that is exponential in the number of entries of the vectors u and i. As we 
have seen above, solving the linear complementary problem plays the role 
of performing a cell identification in our framework. The exponential com- 
plexity is related to the fact that we may have to check membership of all 
cells when evaluating the piecewise linear mapping. However, once a feasi- 
ble set of complementary vectors u and i has been found only one pivoting 
operation is needed in order to determine the new set-up when a single 
constraint has been violated. This has allowed the development of fast and 
memory efficient simulation programs based on this model [38]. 



2.6 Comments and References 

Many important remarks can be made to the developments described so far. 
Rather than obstructing the general presentation with long discussions, we 
have chosen to collect such remarks in a special section at the end of each 
chapter. Some of the material presented in these sections are small remarks, 
while other material discusses related work, gives alternative perspectives 
on the material or presents issues that are not otherwise covered in the book. 

Piecewise Linear or Piecewise Affine? 

The term piecewise linear may at first appear inappropriate for the system 
(2.3) since the dynamics is in fact affine in the state. However, since the name 
is generally accepted, we have chosen not to make a stronger point out of 
this. One may motivate the name piecewise linear by the fact that around 
any trajectory inside a cell, the dynamics will behave linearly. 

Special Classes of Piecewise Linear Systems 

Systems with piecewise linear dynamics have been studied for a long time 
in science and engineering, and many particular classes of piecewise linear 
systems have been introduced. 
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Linear systems with saturation and relay nonlinearities probably have 
the longest tradition (see, e.g., [75, 153, 3]). They are practically relevant and 
serve as excellent prototype systems for understanding more general piece- 
wise linear systems. However, even these apparently simple systems can 
exhibit astonishingly complex behaviours, and the set of analysis and syn- 
thesis tools is still not complete (see [67, 95] for recent developments). 

More general classes of piecewise linear systems, notably the systems de- 
scribed on the implicit form (2.16), have been studied in the circuits and sys- 
tems community [198, 127]. Complementarity conditions also occur in the 
modeling of impulse and contact forces in mechanics leading to dynamical 
systems that are conveniently described by the implicit model (2.16) (see, 
e.g., [131]). Recently, a number of important system-theoretic results have 
been developed for systems on the form (2.16) under the name linear comple- 
mentarity systems [202, 79]. 

Several additional piecewise linear models have been suggested for sys- 
tems with discrete-time dynamics, notably the min-max-plus-scaling sys- 
tems [183] and the mixed logical dynamical systems [15]. A nice exposition 
of the different system models and their relationships can be found in [80]. 

Well-Posedness of Piecewise Linear Systems 

We have defined trajectories of piecewise linear systems as absolutely con- 
tinuous functions that satisfy the model equations at all times. Although any 
reasonable model of a physical system should produce solutions with this 
property, we have seen how model simplifications can result in systems with 
attractive sliding modes. It is thus natural to ask if a given piecewise linear 
system is well-posed in the sense that it produces unique trajectories from 
every initial state. Conditions for well-posedness of bi-modal piecewise lin- 
ear systems (without affine terms) have been derived in [87]. Further results 
have been developed for particular classes of piecewise linear systems such 
as relay systems or linear complementarity systems, see e.g, [95, 200, 78]. As 
we have illustrated by several examples, conditions for well-posedness of 
general piecewise linear systems appears to be very hard to establish. 

Constructing Piecewise Linear System Models from Data 

In one of the motivating examples, we constructed a piecewise linear ap- 
proximation of a smooth nonlinear system by evaluating the right-hand-side 
of the smooth differential equation at a number of points and using linear 
interpolation between these points. This is also how piecewise linear circuit 
models were constructed in [88, 149, 40]. In many cases, however, we do 
not have access to an explicit system model, but need to estimate our model 
from data. Methods for estimating piecewise linear models from data have 
been developed in, for example, [22, 145, 188, 35, 128, 108, 173]. Methods for 
identifying piecewise linear sector bounds from data are described in [107]. 
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3 

Structural Analysis 



In this chapter, we will present some basic tools for analysis of piecewise 
linear dynamical systems. We will treat equilibrium computations and static 
gain analysis as well as verification of affine state constraints and detection 
of attractive sliding modes. These results are useful for ruling out degen- 
eracies in piecewise linear models, provide critical engineering insight, and 
are valuable complements to the Lyapunov-based methods developed in the 
subsequent chapters. 



3.1 Equilibrium Points and the Steady-State Characteristic 

An initial problem in the study of a nonlinear system is to determine its 
equilibrium points. We will understand the term equilibrium point as a con- 
stant trajectory (in the sense of Definition 2.1). Contrary to a linear system, 
which always has an equilibrium at the origin, a general nonlinear system 

x = f(x,u) 
y = g{x , u) 

may have any number of isolated equilibrium points. To compute the equi- 
librium points, we pick a constant input u = u eq and find the solutions x eq 
to the equation f(x eq , u eq ) = 0. By repeating the computations for a range of 
inputs and recording the stationary outputs y eq = g(x eq , w eq ) we obtain the 
steady-state characteristic . 

Equilibrium Point Computations 

Computing the equilibrium points of a piecewise linear system is straight- 
forward: we simply have to check if each affine subsystem has any equilib- 
rium points within its operating regime. Hence, for given u eq , we need to 
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3.2 Equilibrium Points and the. Steady-State Characteristic 
chock, if there exists any solutions :r eq e to 

0 — -A^X^q + &i H" 

The following result shows how the equilibrium point computations can be 
conveniently formulated as a linear programming problem. 

Proposition 3.1— Equilibrium points 

The set of equilibrium points of the piecewise linear system (2.3) in the inte- 
rior of cell A 'i is given by 

Q.i = {x | AiX + in + Bi'Uoq = 0. Gix + cn >- 0 j 

Finding a feasible equilibrium point :r eq € Qi, or verifying that no such 
point exists, is a linear programming problem. □ 

The computations can be simplified when the matrices A* are invertible. In 
this case the potential equilibrium point x oq = —A { 1 (<ii + B}U 0 q ) lies in the 
interior of X x if and only if 

— GiA^ l (Oi | I $i>- 0. 

The above approach is still computationally intensive since all cells have 
to be considered one-by-one. More efficient algorithms have been developed 
for piecewise linear systems with continuous right-hand side, see [44, 207, 
155]. The following example demonstrates the analysis. 

Example 3.1— A Bistable Electrical Circuit 

Consider the tunnel diode circuit in Figure 3.1. The circuit equations are 



c 5T- 


r. f-L ~ Qr(Vc.) 


T <l ■ 

l m“ 


, = u - Ri L - % 



i>L L 




Figure 3.1 Bistable electrical circuit (left) 2 




piecewise linear characteristic (right). 
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Here, g R (v c ) denotes the conductance of the nonlinear element. We use the 
parameters R = l.Skfi, C = 2pF, L = 5 fiH (taken from [116, 43]). By mea- 
suring time in nanoseconds and currents in mA, we obtain 

= 0.5 [-g R (y c ) + ^l] 

^-i L = 02[-v c - Ri L + u] 

A piecewise linear approximation of the system is constructed by approxi- 
mating g R (v c ) by a seven-segment piecewise linear function, see Figure 3.1 
(right). Letting u eq =1.2, we run the computations of Proposition 3.1 to find 
the three equilibrium points 



x 



* 

l 



0.07 
0.76 ’ 



x 2 = 



0.28 
0.62 ’ 



x 3 = 



0.87 

0.22 



which have good correspondence with the computations in [116]. □ 

The value of an equilibrium point computation is significantly increased if it 
is combined with a local analysis of the properties of the equilibria. The lo- 
cal analysis, which can be done by eigenvalue inspection of the matrices A it 
might reveal important qualitative properties of the overall system. For ex- 
ample, if the system has no sliding modes and no stable equilibrium points 
then the state vector either tends to infinity or towards some non-stationary 
behavior (possibly a limit cycle). 

The Steady-State Characteristic 

The extension to computation of the steady-state characteristic is now im- 
mediate. For each admissible input u eq , we simply compute the associated 
equilibrium points x eq (as we have seen above, there might be many) and 
record the corresponding outputs 

Ueq = C^Xq q T~ Ci T~ DiU eq 

This results in a steady-state characteristic {w eq , {?/eq}}- 

Example 3.2— Steady-State Characteristic of Bistable Circuit 
As shown in Example 3.1, the bistable electrical circuit exhibits multiple 
equilibrium points for u eq = 1.2V. By considering the capacitor voltage v c as 
an output and performing the static analysis above with u eq e [0, 2], we ob- 
tain the characteristic shown in Figure 3.2. We see that the system has three 
equilibrium points for 0.6V < u eq < 1.6V, while it has unique equilibria for 
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Figure 3.2 Steady-state characteristic showing multiple equilibria for 0.6V < u < 
1.6V. Local stability of the equilibria is marked by o, instability marked by x. 



other inputs. One equilibrium point is locally unstable (marked x) while the 
others are locally stable (marked o). Indeed, the circuit has been used as a 
computer memory, exploiting the fact that a change in the input voltage can 
shift the system state from one equilibrium to the other. □ 



State Transformations 

When analyzing the global properties of an equilibrium point, it is often 
convenient to make a state transformation so that this point is the origin in 
the transformed coordinates. This transformation is then given by 
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x := Tx. 



By applying the state transformation x \ —> T 1 z, one obtains a piecewise 
linear system with an equilibrium point at the origin. 



3.2 Constraint Verification and Invariance 

In practice, it is often desirable to verify that a system satisfies safety con- 
straints in the presence of exogenous disturbances. In this section we will 
show how a variation of this problem can be studied using vector field con- 
siderations and resolved using linear programming. 

To this end, consider the piecewise linear system (2.3) and assume that 
u(t) is an exogenous disturbance that satisfies the magnitude constraint 

u<u(t)<u forall£>0 (3.1) 
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but is otherwise arbitrary. Let the set of safe operation be given by 

O = {x | Rx + s y 0} = {x | r^x + s k > 0 k = 1 , . . . , m} (3.2) 

We are interested in verifying that all trajectories that originate in O remain 
in the set for all future times. We then say that the set O is positively invariant 
with respect to the dynamics (2.3) (see, e.g., [116, 24]). 




Figure 3.3 Desired invariant set O (left). Invariance of the individual constraints 
r^x + s k > 0 requires that the vector field is inward at the constraint hyperplane. 

Invariance of O can be addressed by establishing invariance of each in- 
dividual constraint separately. For each constraint 

r[x + Sk >0 (3.3) 

in (3.2), define the associated set of constrained states Hk as 



Hk = {x £ O | rjx + s k = 0} 



Then, the constraint (3.3) is invariant if and only if 

^( r k x (t) + Sk) > 0 for all x{t)eHk (3.4) 

This condition ensures that if a trajectory approaches the constraint set from 
within O, the quantity rjx(t) + s k will increase in value, and the constraint 
(3.3) will remain being satisfied. 

To express the invariance condition (3.4) for a piecewise linear system, 
define I k as the set of indices for the cells that intersect the constraint set, 

I k = {i | Xi fl Hk 7^ 0}. 
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3.3 Detecting Attractive Sliding Modes on Cell Boundaries 

For ease-of-computation, we will require the inequality in (3.4) to be strict. 
Then, to verify (3.4), we need to verify that 

r^(AiX + di + Biu) > 0 for all x € Xi fl H k , i € h, (3.5) 

and all u e [u,u\. The geometrical interpretation is that the vector field 
should be inward at Hk, see Figure 3.3. Condition (3.5) is equivalent to es- 
tablishing that there is no solution x to the inequalities 
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A direct application of Farkas’ lemma now yields the following result. 
Proposition 3.2— Constraint Verification 

Consider the system (2.3) whose input satisfies the constraint (3.1). If there 
exists vectors t >z 0 and w iki h 0 that satisfy 

v ik^ik = 0 v ik s ik{lk) > 0 

w Ik S ik = 0 wf k S ik (u) > 0 

for every i e I s and k e {1, . . . , m}, then O is positively invariant. □ 

For related results on invariance in piecewise linear systems, see e.g., [39, 49, 
92, 173] and the Comments and References section at the end of this chapter. 



3.3 Detecting Attractive Sliding Modes on Cell Boundaries 

As discussed in Chapter 2, piecewise linear systems with discontinuous 
right-hand side may exhibit sliding modes on the boundaries between two 
or more cells. Since the presence of sliding modes complicates the system 
analysis and challenges the assumptions made in the modeling process, it is 
important to have tools that detect the presence of attractive sliding modes 
in a given model. A complete, yet computationally efficient, solution to this 
problem appears to be out of reach. We have seen how sliding modes can 
occur on the intersection of multiple boundaries even if the vector fields are 
not traversal, in which case one needs to do a full stability analysis to assess 
if the sliding mode is attractive. We will therefore limit our ambition and 
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only consider attractive sliding modes on the n — 1 dimensional boundaries 
between two cells. Initially we consider the case of a regular sliding mode, 
where the vector fields in both cells point toward the common boundary as 
in Figure 3.4. The following definition is then natural. 

Definition 3.1— Regular Sliding Set 

Let dXij = Xi n Xj C {x | hJjX + gij = 0} be the boundary between cells 
Xi and Xj, defined such that hfjX + gij > 0 for x G X The set 

Sij = {x e int dX^ | hfj(AiX + a*) < 0 A hfj(AjX + dj) > 0} 

is called the regular sliding set for (2.3) on dX ^ . If the set is non-empty, then 
the system (2.3) has an attractive sliding mode on S^. □ 

The following proposition, which shows how attractive sliding modes can 
be detected by solving a simple linear program, is now immediate. 

Proposition 3.3 

If the optimal solution (x*, t*) to the linear programming problem 
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has t* > 0, then (2.3) has a non-empty regular sliding set on dX^. □ 

A regular sliding mode occurs when the vector fields in both cells point 
toward the common boundary dX^ = {x e Xi fl Xj \ hJjX X g%j = 0}. Slid- 
ing modes can also arise when the vector fields are tangent to the bound- 
ary, which gives rise to so called higher-order sliding modes [60]. For example, 
a second-order sliding mode occurs when the first time-derivative of hJjX 
vanishes on some point on the boundary 

hJj(AiX + Oi) = hfj(AjX + <ij) = 0 for some x G dX^ 

while the second derivative of hj-x is negative for the dynamics in Xi and 
positive for the dynamics in Xj. A third-order sliding mode occur when the 
first two derivatives of hJjX vanish, etc. Similarly to above, these conditions 
can be used to define higher-order sliding sets and the existence of such a 
sets on cell boundaries can be verified using linear programming. The re- 
sulting conditions are closely related to the well-posedness test for bimodal 
piecewise linear systems proposed in [87]. 



38 




3.4 Comments and References 




dXij 




Figure 3.4 Vector fields in neighboring cells point towards the regular sliding set Sij. 



3.4 Comments and References 



From Signals to Symbols 

Vector field considerations are useful for obtaining simplified pictures of the 
dynamics of piecewise linear systems. By examining the behavior of trajec- 
tories on the surfaces of the cells it is sometimes possible to limit the sys- 
tem behavior to a finite number of alternatives. For example, in some cases 
one may be able to establish that no trajectory that start inside a cell can 
exit through a certain cell face. Such an analysis gives a natural aggregation 
of solutions that makes it possible to abstract away detailed information 
about solutions in order to obtain a simple picture of the global dynamics. 
This idea has been used in [159] for the construction of “phase-portraits” for 
high-dimensional piecewise linear systems. Similar ideas have been used in 
[92] for verification of piecewise linear hybrid systems. Related is also the 
concept of cell-to-cell mappings, see [84]. 

Controllability and Observability of Piecewise Linear Systems 

The treatment of static gain analysis touches upon the concepts of observ- 
ability and controllability of piecewise linear systems. These issues have not 
been investigated within this book. Controllability of a certain class of piece- 
wise linear systems has been treated in [208, 125]. Related is also the more 
recent work [203]. 

Constraints and Invariance 

Invariant sets is a useful notion in the analysis and design of control systems. 
Their history in control dates back to the early work of Lyapunov: as any 
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level set of a proper Lyapunov function is invariant, Lyapunov functions can 
be used to establish invariance of state and control constraints [24, 117, 27]. 
More recently, invariance have been used in the design of controllers that 
minimize the peak-to-peak gain -induced gain) see [25, 186] and the ref- 
erences therein. These approaches are often based on the same ideas that 
were used in Section 3.2. Note, however, that the problem of computing an 
invariant set is much more involved than simply verifying that a given poly- 
hedral set is invariant. An efficient approach for computing invariant sets of 
piecewise linear systems has been developed in [92]. 
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4 

Lyapunov Stability 



The main contribution of this book is the development of a Lyapunov-based 
analysis method for piecewise linear systems. The key component of such 
an analysis, namely methods for Lyapunov function construction, will be 
presented in this chapter. We will show how piecewise quadratic and piece- 
wise linear Lyapunov functions can be computed via convex optimization. 
These approaches are substantially more powerful than analysis based on 
quadratic Lyapunov functions, yet allow the analysis to be be carried out 
using efficient numerical computations. The application to system analysis 
and design of optimal control laws is described in the subsequent chapters. 

As always, it is a good idea to use “simple things first”. After determin- 
ing the equilibrium points of a piecewise linear system it is advisable to ver- 
ify local stability properties first. For locally stable equilibrium points, one 
may invoke the tools developed in this chapter to try to extend the domain 
of analysis and in many cases establish global results. 



4.1 Exponential Stability 

Stability is one of the most fundamental properties of dynamical systems. 
Intuitively, stability is the property that a system does not “blow up”. In 
many control systems, stability alone is not satisfactory but one seeks to 
achieve asymptotic stability, which ensures that the system state comes to 
rest after an initial transient. We will focus on the particular case of expo- 
nential stability, which guarantees that the convergence of the system state 
to its equilibrium point can be bounded by an exponential function of time 
(see, e.g., [116] for precise definitions). 

There are certainly many asymptotically stable systems whose conver- 
gence is not exponential. Still the framework of exponential stability is par- 
ticularly attractive for Lyapunov-based analysis of piecewise linear systems. 
For smooth nonlinear systems, an equilibrium point is locally exponentially 
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stable if and only if its linearization around this point is exponentially sta- 
ble ([116, Theorem 3.13]). In other words, the linearization provides neces- 
sary and sufficient information about local exponential stability. Moreover, a 
smooth nonlinear system is globally exponentially stable if and only if there 
exists a Lyapunov function that proves it ([116, Theorem 3.12]). 

The following result, which combines a number of standard results from 
Lyapunov theory, will be the main tool throughout this chapter. 

Lemma 4.1 

Let x(t) : [0, oo) — » M n and let and let V(t) : [0, oc) — » R be a non- 
increasing and piecewise C 1 function satisfying 

J t v (t) < -^\x(t)\\ p (4.1) 

for some 7 > 0 and some p > 0, almost everywhere on [0, oc). 

If there exists a > 0 such that 

a\\x(t)\\ p <V(t)< P\\x(t)\\ p (4.2) 

then \\x(t)\\ tends to zero exponentially. If the maximal a that satisfies (4.2) is 
negative, then ||x(£)|| — ► 00 as t — ► oc. □ 

Proof: See Section B.2. 

Note that the above result allows both verification of exponential stabil- 
ity and detection of instabilities. Moreover, the formulation does not require 
that V(t) be continuous as long as the value of V(t) decreases at the points 
of discontinuity. 



4.2 Quadratic Stability 

The term quadratic stability refers to stability that can be established using 
a quadratic Lyapunov function. Quadratic stability dates back to the pio- 
neering work of Lyapunov [133, 134] who established that the existence of 
a quadratic Lyapunov function is a necessary and sufficient condition for 
asymptotic stability of a linear system. Quadratic Lyapunov functions are 
often the first resort also in the analysis of nonlinear systems and much 
work on absolute stability is based on quadratic Lyapunov functions (see, 
e.g., [45, 27]). 
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Quadratic Stability for Local Linear Analysis 

To verify exponential stability of a nonlinear system one may start by es- 
tablishing stability of its linearization. This is often done by eigenvalue in- 
spection. One drawback of this approach is that one has no idea about the 
domain of validity of the analysis and exponential convergence can only be 
guaranteed for trajectories starting very close to the equilibrium. A useful 
addition to a local stability analysis is to also have an estimate of the domain 
of attraction of the equilibrium point, i.e., the set of initial values from which 
trajectories are guaranteed to converge to the equilibrium. In the piecewise 
linear framework the domain of validity of the linearization is incorporated 
in the model. By using Lyapunov function computations rather than eigen- 
value inspection it is then possible to compute a guaranteed domain of at- 
traction. The following theorem establishes exponential stability of an equi- 
librium point and also returns the largest domain of attraction that can be 
guaranteed by a local analysis using quadratic Lyapunov functions. 

Proposition 4.1— [27] 

Let x = Ax be valid in the polyhedron X 0 = {x \ g^x < 1 7 k = 1, . . . ,p}. 
The origin is exponentially stable if and only if the convex optimization 
problem 



minimize log det Q 1 
subject to 0 < Q 

0 > QAA \- AQ 

1 > 9 kQ T 9 k k = l ,..., p . 

has a solution. Moreover, the ellipsoid £ roa = {x \ x T Q~ 1 x < 1} is the do- 
main of attraction with largest volume in X 0 that can be estimated using any 
quadratic Lyapunov function. □ 

Proof: See [27, Section 5.2.2]. 

The strength of this result is that the computations, in addition to ver- 
ifying exponential stability, return the largest region of attraction that can 
be estimated using any quadratic Lyapunov function. In most cases, this do- 
main is substantially larger than what may be obtained by simply solving 
the Lyapunov equation 



A t P + PA = -R (4.3) 

for some arbitrary chosen positive definite matrix R. This is illustrated in 
the following example. 
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Example 4.1— Local Analysis of Saturated System 
Consider the saturated linear system of Example 2.1, 



and let 



x = Ax + bsat(k T x) 



0 


l" 


, b = 


"o" 


, k = 


-2 


0 


0 




1 




-3 



By restricting the analysis to the unsaturated region, the polyhedron \k T x\ < 
1, we can use Proposition 4.1 to verify exponential stability of the origin. 
The analysis returns the region of attraction shown in dashed lines in Fig- 
ure 4.1. If we simply solve the Lyapunov equation (4.3) with R = /, we can 
only guarantee exponential stability of trajectories starting within the (much 
smaller) domain shown as the filled ellipsoid in Figure 4.1. □ 




Figure 4.1 Linear (white) and saturated (shaded) regions in the state space. The larger 
region of attraction is computed via Proposition 4.1, while the smaller (filled) ellipsoid 
is the stability domain obtained by solving the Lyapunov equation (4.3) with R = I. 

Although Proposition 4.1 guarantees a much larger region of attraction than 
the Lyapunov equation approach, the result is still very weak. In fact, the sat- 
urated system is globally asymptotically stable (we will prove this in Chap- 
ter 7, where we also develop specialized methods for estimating regions 
of attraction for linear systems with saturation). One of the reasons for the 
disappointing performance is that we have only analyzed the behaviour in 
the unsaturated region. It is possible that the computed Lyapunov function 
would be able to verify stability also for the saturated regions, but since 
these have not been accounted for in the analysis no such conclusion can be 
drawn. This situation will appear repeatedly when we compute Lyapunov 
functions on a restricted domain: stability can only be granted for trajec- 
tories whose initial values lie within the largest level set of the computed 
Lyapunov function that is fully contained in the analysis region. 
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Quadratic Stability for Global Nonlinear Analysis 

Quadratic Lyapunov functions are often the first resort also for analysis of 
nonlinear systems. Many of these methods are somehow related to the anal- 
ysis of the linear differential inclusion 



x e co{Aix, . . . ,Alx}. (4.4) 

In other words, they are related to the analysis of the family of linear time- 
varying systems that can be written as 



L 

x = A(t)x = oti{t)Aix{t) (4.5) 

i= 1 



where ai(t) are arbitrarily time-varying weights that satisfy ai(t) > 0 and 
Yli=i a i(f) = 1- The following result is central [83, 30]. 



Proposition 4.2 

Consider the system (4.4). If the convex optimization problem 



0 <P = P T 

0 > Aj 1 P T- P Ai i = 



(4.6) 



has a solution, then the origin is globally exponentially stable. □ 

Proof: See Section B.2. 



Note that Proposition 4.2 only gives sufficient conditions for stability. 
Since the Lyapunov function search in Proposition 4.2 is a convex optimiza- 
tion problem a solution P can always be found if it exists. The conservatism 
comes from the fact that quadratic Lyapunov functions are only sufficient 
for establishing stability of linear differential inclusions [142]. 

The main advantage of Proposition 4.2 is that the search for a quadratic 
Lyapunov function has been formulated as a convex optimization problem. 
The computational complexity for solving the inequalities in Proposition 4.2 
grows gracefully with the number of extreme systems L, and solving the 
multiple Lyapunov inequalities is not much more demanding than solving 
a single Lyapunov equation, see [205, 27] 

There have been many applications of the above results to systems with 
piecewise linear dynamics, see for example the work on fuzzy systems [195, 
219] which can be embedded in the linear time varying formulation (4.5) 
by the appropriate restrictions. The application to piecewise linear systems 
typically takes the following form. 
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Corollary 4.1 

Consider the system (2.3), and assume that ai = 0 for every i G /. If the con- 
vex optimization problem (4.6) has a solution, then every trajectory x(t) e 
U i e iXi of (2.3) with u = 0 tends to zero exponentially. □ 

Proof: See Section B.2. 



In some cases, it is of interest to verify that no common solution P to the 
conditions of Proposition 4.2 exists. This verification can be done by solving 
the following dual problem (cf [14]). 



Proposition 4.3 

If there exists positive definite matrices Ri satisfying 

{ 0 < Ri = Rf 

L for i ^ 

0 < ^^( RjAj yAiRi) 

i= 1 

then there exists no solution to the LMIs of Proposition 4.2. 
Proof: See Section B.2. 



(4.7) 

□ 



4.3 Conservatism of Quadratic Stability 



Although computationally attractive, several issues make the quadratic sta- 
bility analysis of piecewise linear systems given in Corollary 4 . 1 very conser- 
vative. A first limitation is that no affine terms are permitted in the dynamics 
so that a simple system such as the saturated control system in Example 4.1 
can not be analyzed as it stands. A second issue is that no information about 
the partition is used in the analysis. Rather than exploiting that each sys- 
tem of differential equations only describes the system dynamics in a re- 
stricted part of the state space, the piecewise linear dynamics is embedded 
in a global differential inclusion. The following example illustrates the limi- 
tations of this approach. 



Example 4.2— The Need for Using Partition Information 
Consider the piecewise linear system 



/ 


’- 0.1 1 




-10 - 0.1 


< 






’- 0.1 10 " 


< 


-1 - 0.1 



X\X2 > 0 



X\X2 < 0 



(4.8) 
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The system matrices are stable and have the same eigenvalues. The simula- 
tion shown in Figure 4.2 indicates that the system is stable. 

It is straightforward to verify that V (x) = x T x is a Lyapunov function 
for the system. Still, an application of Proposition 4.2 fails. Since the matrices 



Ri = 



2 

-1 



-1 

2 



R2 



2 1 
1 2 



satisfy the conditions in Proposition 4.2 this failure is fundamental (and not 
an artifact of numerical optimization) . We can understand the reason for this 
failure by interchanging the system matrices in the model (4.8). Simulating 
this system yields unbounded trajectories, see Figure 4.2. Since stability de- 
pends on the partition, Proposition 4.2 cannot be used to establish stability. 
At the end of this chapter we will return to this example with a more pow- 
erful tool set and prove stability of the initial setup. □ 










) 



o iio 



Figure 4.2 The original setup (left) is exponentially stable, while interchanging the 
two system matrices gives an unstable system (right; note the different scalings!) . Thus, 
stability can not be proved without taking the structural information into account. 



A third limitation is of course that many systems do not admit a quadratic 
Lyapunov function. The following example illustrates a simple system that 
cannot be analyzed using quadratic Lyapunov functions. 

Example 4.3— The Need for Non-Quadratic Lyapunov Functions 
C onsider the piecewise linear system x(t) = Aix(t) with the cell partition 
shown in Figure 4.3 and system matrices 



II 

CO 

II 


—e co 


to 

II 

II 


— e aco 




—aco —e 




—co —e 
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Letting a = 5, co = 1 and e = 0.1, a simulation with initial value xo = (—2, 0) 
generates the trajectory shown in Figure 4.3. Clearly, no quadratic Lyapunov 
function can generate level sets with the property that once a trajectory en- 
ters a level set, it remains in this set for all future times. Hence, there is no 
obvious quadratic Lyapunov function that can establish asymptotic stability. 

-^2 




Figure 4.3 Simulated trajectory (full) and cell boundaries (dashed) in Example 4.3. 



□ 

As illustrated above, the computational convenience of quadratic stability 
analysis of piecewise linear systems is out-weighted by its many limitations. 
The analysis conditions do not permit affine terms in the dynamics, and 
even very basic piecewise linear systems cannot be treated as they stand. 
The use of partition information can be critical for establishing stability, yet 
it is completely disregarded in Corollary 4.1. Finally, it can be very conser- 
vative to only consider quadratic Lyapunov functions. 

In the remaining parts of this chapter we will develop an approach that 
does not suffer from these shortcomings. The method will be based on non- 
quadratic Lyapunov functions, exploit partition information and allow for 
affine terms in the dynamics. All analysis conditions will be formulated as 
convex optimization problems, allowing complex piecewise linear systems 
to be analyzed using efficient numerical computations. 



4.4 From Quadratic to Piecewise Quadratic 

To find inspiration for alternatives to the globally quadratic Lyapunov func- 
tions, we will analyze the simple selector control system shown in Figure 4.4. 
Selector control is a common strategy for constraint handling in the process 
industry (see, e.g, [8, 59, 31]) that often results in closed-loop systems with 
piecewise linear dynamics. 

Consider the set-up shown in Figure 4.4 (left) and assume that Go is a 
linear system described by the state-space equations x = Aqx + Bu. The 
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Figure 4.4 Selector control system (left) transformed into feedback form (right). 



closed-loop dynamics can then be written as 

x = Aqx + B max (kjx, k^x) = (Ao + Bk^)x + B max ((Aq — k2) T x : 0) 
Letting A = A 0 + Bk J and k = ki - fc 2 , we can re-write the dynamics as 

x = Ax + B max (k T x, 0) (4.9) 

The system is clearly piecewise linear, and can be written as 

( Ax if k T x < 0 

\(A-\-Bk T )x if k T x > 0 
Now, consider the specific system defined by 



"-5 -4 


, B = 


V 


, km 


Y 


-1 -2 




21 




0 



By solving the dual problem stated in Proposition 4.3 one can verify that 
there is no globally quadratic Lyapunov function V (x) = x T Px that verifies 
stability of the system. Still, the simulations shown in Figure 4.5 indicate 
that the system is stable. 

As an alternative to a globally quadratic Lyapunov function it is natural 
to consider the following Lyapunov function candidate 



V(x) 



x T Px, if k T x < 0 

x T Px + r}(k T x) 2 , ifk T x>0 



(4.10) 



where P and rj € R are chosen so that both quadratic forms are positive def- 
inite. Note that the Lyapunov function candidate is constructed to be con- 
tinuous and piecewise quadratic. The search for appropriate values of r\ and 



49 










Chapter 4. Lyapunov Stability 





X ‘} 4 



Figure 4.5 Trajectories in the phase plane of the selector control system. 



P can be done by solving the following linear matrix inequalities 

P = P T >0, P + rjkk T > 0, 

A t P + PA <0 (A + Bk T ) T (P + rjkk T ) + (P + (^4 + T?/c t ) < 0. 

One feasible solution is given by P = diag{l, 3} and rj = 9, which estab- 
lishes exponential stability of the system. The level surfaces of the computed 
Lyapunov function are indicated in Figure 4.5. 

Relation to Frequency Domain Criteria 

It is instructive to compare this solution with what can be achieved using 
frequency domain methods such as the circle and Popov criteria. In order to 
put the system in negative feedback form, we re-write (4.9) as 

x = Ax - (-B) max (k T x, 0) (4. 1 1) 

Defining G(s) = —k T (si - A)~ 1 B, we obtain the frequency condition 

Re G(iuj) > —1, Vte G [0, oo] 

for the circle criterion and 

Re [(1 + icor])G(iu;)\ > —1, Vu; G [0, oo] 

for the Popov criterion. Inspection of the Nyquist and Popov plots of Fig- 
ure 4.6 reveals that stability follows from the Popov criterion but not from 
the circle criterion. 

The failure of the circle criterion comes as no surprise, since it relies on 
the existence of a common quadratic Lyapunov function (see, e.g., [116]) 
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Figure 4.6 The circle criterion (left) fails to prove stability. The Popov plot (right) is 
separated from —1 by a straight line of slope 1/77. Hence stability follows. 



which we know does not exist. The Popov criteria, on the other hand, is 
based on the Lyapunov function 



V(x) 



x T Px -\~2rf y?(cr) da. 

Jo 



(4.12) 



By evaluating this function for <p(a) = max (a, 0) and z = k T x, one recovers 
the Lyapunov function candidate (4.10) used in the optimization above. 

It is not hard to establish that all Lyapunov function of Lure-type 



V(x)=x T Px-\- f <p(a) da 

Jo 

constructed from piecewise linear functions <p(a) are continuous and piece- 
wise quadratic. However, rather than tailoring the analysis to linear sys- 
tems interconnected with scalar nonlinearities we will aim for results that 
are applicable to general piecewise linear systems. Motivated by the exam- 
ples above, we will employ Lyapunov functions that are continuous and 
piecewise quadratic. 



4.5 Interlude: Describing Partition Properties 

In order to generalize the ideas used in the motivating example we need 
to be able to construct continuous piecewise quadratic functions on gen- 
eral polyhedral partitions and we need to be able to enforce positivity (and 
negativity) of such functions. In this section, we will introduce a matrix pa- 
rameterization of continuous piecewise quadratic functions and describe a 
computational procedure for enforcing positivity. 



51 



Chapter 4. Lyapunov Stability 

A Parameterization of Continuous Piecewise Quadratic Functions 

Our analysis will be based on Lyapunov functions that are continuous and 
piecewise quadratic. A fundamental problem in the search for such func- 
tions is how the continuity constraint should be enforced on the Lyapunov 
function candidate. One potential solution would be to patch together piece- 
wise quadratics “by hand”, as was done in the analysis of the selector sys- 
tem. If we only require continuity of the piecewise quadratics, we can allow 
less restrictive updates in the quadratic forms than pkk T used above. To il- 
lustrate this, let 



V (x, i) = x T PiX + 2 qj Xi + n 



X 


T 


Pi Qi 




X 


1 




1 




1 



X T PiX 



be the quadratic expression used to describe V(x) for x e Xi. Consider two 
cells Xi and Xj with boundary hyperplane dHij = {x \ hLx + gij = 0}. To 
obtain continuity across the boundary hyperplane 



V(xJ) = V(x,i) Vx G dHij 

there must exist Lj e R n and £ R such that 



V(x,j) = V(x,i)+2(hf j x + g ij )(tT J x + s ij ). 

Letting h+j = [hj- gij\ T and Lj = \tj- s^] T this condition reads 

Pj = Pi + + iij hjj . (4.13) 

Since each cell boundary induces a new equality constraint on the form 
(4.13) it is very cumbersome to construct explicit expressions for the matri- 
ces Pi other than in very simple examples. An alternative is to simply view 
(4.13) as an equality constraint in the variables P it Pj and Lj that the op- 
timization software should deal with. However, as of today, very few SDP 
solvers treat equality constraints (see the Comments and References section 
at the end of this chapter for more details). 

We will therefore introduce a compact matrix parameterization of contin- 
uous piecewise quadratic functions on polyhedral partitions. This will allow 
our analysis procedures to be implemented in widely available SDP solvers. 
The parameterization is based on continuity matrices, defined as follows: 

Definition 4.1— Continuity Matrix 
A matrix Fi = [Fi fi] is a continuity matrix for cell Xi if 

Fix(t) = Fjx(t) for x(t) £Xi n Xj. (4.14) 
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Furthermore, we say that {Fi} have the zero interpolation property if 

fi = 0 forie/o. (4.15) 

□ 

A compact matrix parameterization of continuous piecewise quadratic func- 
tions can now be obtained as follows. 

Lemma 4.2— PwQ Parameterization 

Let {Xi} ie i be a polyhedral partition, and let Fi £ R^ x ( n+1 ) be continuity 
matrices that satisfy (4.14). Then, for each T £ R pxp t the scalar function 

V(x) = x T FjTFiX := x T PiX for x £ X*. (4.16) 

is continuous and piecewise quadratic. Moreover, if {Fi} have the zero in- 
terpolation property, then there exist a and (3 such that 

a\\x\l<V{x)<p\\x\l (4.17) 

□ 



Proof: See Section B.2. 



Note that the continuity matrices for a given partition are not unique. 
For example, if we use the following continuity matrix in all regions 



F = 



J nx 1 



ieL 



the parameterization suggested in Lemma 4.2 results in a Lyapunov func- 
tion candidate that is globally quadratic. Clearly, one would like a way of 
constructing the continuity matrices that gives maximal freedom in the Lya- 
punov function search. Procedures for computing continuity matrices for a 
given piecewise linear system are described in Appendix A. The following 
example illustrates one choice of continuity matrices for the saturated sys- 
tem from Example 2.1. 



Example 4.4— Continuity Matrices for Saturated System 

The following matrices are natural continuity matrices for the saturated 

feedback system in Example 2.1. 





'- k T 


-1 




Olxn 


o" 




Olxn 


0 " 


Fi = 


Olxn 


0 


II 


Olxn 


0 


, A = 


k T 


-1 




X 

to 


0 




hx2 


0 




hx2 


0 
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Note that the matrices have the zero interpolation property and that a Lya- 
punov function candidate constructed as in Lemma 4.2 has no affine terms 
in the region that contains the origin. □ 



Verifying Positivity of a Quadratic Function on a Polyhedron 

In order for a function to satisfy the conditions of Lemma 4.1, we need to 
enforce that the function is positive and that its time-derivative is negative. 
For piecewise quadratic functions, the positivity condition reads 

x T PiX > c*||x ||2 Vx G Xi , i G I (4.18) 

Hence, for each i we need to verify that there exists a positive scalar a such 
that x T PiX— ax T x is positive on Xi. We know from Example 4.2 that it is crit- 
ical to exploit the partition information and that it is too restrictive to require 
that the expression is positive for all x G M n . To support the computations, 
we define polyhedral cell boundings as follows. 

Definition 4.2— Polyhedral Cell Bounding 
A matrix Ei = [Ei e*] is called a polyhedral cell bounding if it satisfies 

Eix(t) y 0 forx(t) G Xi. (4.19) 

Furthermore, we say that {Ei} have the zero interpolation property if 

ei = 0 for i G Jo. (4.20) 

□ 

Verifying positivity of a piecewise quadratic function on a polyhedral parti- 
tion can now be done using the following result. 

Lemma 4.3 

Consider the piecewise quadratic function 

{ x T PiX x G X^ i G Jo 

rp ~ (4-21) 

x PiX x G Xi, i G 1 1 

with Pi = P- T , Pi = Pj, and let Ei be cell boundings satisfying (4.19) and 
(4.20). If there exist matrices W* with nonnegative entries such that 

Pi — Ej 1 WiEi >0 i G Jo 

Pi ~ Ej WiEi > 0 
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i G I\ 



(4.22) 

(4.23) 




4. 6 Piecewise Quadratic Lyapunov Functions 



then there exists a > 0 such that V (x) > a\\x\\\ iov &\\ x e X . □ 

Proof: See Section B.2. 

The technique used in Lemma 4.3 is known as the 5-procedure and is 
treated in more detail in Appendix A. The method is conservative in general, 
but appears to work very well in practice. 

There is a direct relationship between cell identifiers and cell boundings. 
The key difference is that cell boundings use linear forms for bounding cells 
that contain the origin. This will be critical for formulating analysis pro- 
cedures using strict LMIs. In Appendix A we give a simple algorithm for 
translating cell identifiers to cell boundings and show that this step does not 
introduce any conservatism in the Lyapunov function search. The following 
example illustrates cell boundings for the saturated system. 

Example 4.5— Cell Boundings for Saturated System 
The following matrices are cell boundings for the saturated system 



k T 


-l" 




0 1 X 71 


o" 




' k T 


-l" 


blxn 


1 


j = 


Olxn 


0 


, E 3 = 


Ol xn 


1 



Moreover, these cell boundings have the zero interpolation property. □ 



4.6 Piecewise Quadratic Lyapunov Functions 

The previous section has laid the grounds for analysis of piecewise lin- 
ear systems using piecewise quadratic Lyapunov functions. In Lemma 4.2 
we proposed a matrix parameterization of continuous piecewise quadratic 
functions and in n Lemma 4.3 we showed how it is possible to express the 
condition that the function be positive on the partition using linear matrix 
inequalities. We will now combine these results to formulate the search for 
piecewise quadratic Lyapunov functions on the form 



f x T PiX 



for x e X h i e Io 



V(x) = 



r- -1 


T 


<- -1 


X 


p. 


X 


1 


A 


1 



= x T PiX + 2 qj X + Vi for X e X h i e I\. 



(4.24) 



We have the following result. 
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Theorem 4.1— Piecewise Quadratic Stability 

Consider symmetric matrices T, Ui and Wi such that Ui and Wi have non- 
negative entries, while 

Pi = F?TFi, 

Pi = FfTFi , 

satisfy 

0 > Af Pi + PiAi + Ef UiEi 
0 < Pi - EjWiEi 

0 > AJ" Pi -\- PiAi + Ej UiEi 
0 < Pi - EjWiEi 

Then every trajectory x(t) G U keiXi satisfying (2.3) with u = 0 for all t > 0 
tends to zero exponentially. □ 

Proof: See Section B.2. 



i e Io 
i e h 

i e Io 
i e h 



In the absence of attractive sliding modes, the above conditions ensure 
that (4.24) is a Lyapunov function for the system. If the partition X is pos- 
itively invariant then every trajectory that starts in X tends to zero expo- 
nentially. In particular, if the partition covers the whole state space then the 
system is globally exponentially stable. Even if invariance of X cannot be 
established, any level set of V(x) that is fully contained in the partition is a 
region of attraction for the equilibrium point x = 0. 

With the piecewise quadratic stability theorem at hand we can now re- 
turn to the motivating examples where the standard LMI conditions for 
quadratic stability fail. 



Example 4.6— Piecewise Quadratic where Quadratic Fails - I 
Consider the system (4.8) in Example 4.2. Let Ei denote the cell bounding 
used in quadrant i. Letting 



E\ = —Es 



1 

0 



0 



E2 = —E4 



-1 0 
0 1 7 



and Fi = 





we invoke Theorem 4.1 and find a feasible solution 



V(x) = x T x 



x G Xi, i G I . 
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Hence, stability can indeed be established using a quadratic Lyapunov func- 
tion but it is critical to use partition information in the Lyapunov function 
search. The level curves of the computed Lyapunov function are shown to- 
gether with a simulation in Figure 4.7 (left). □ 





Figure 4.7 Level surfaces (dashed) for the systems in Example 4.2 and Example 4.3 
computed using Theorem 4.1. In both cases, the standard conditions for quadratic sta- 
bility fail while Theorem 4 . 1 verifies stability. 



Example 4.7— Piecewise Quadratic where Quadratic Fails - II 
As a second example, consider the system with flower-like trajectories used 
in Example 4.3. Similarly as above, we let 



Ei = —E3 = 



-1 1 
-1 -1 



E2 = —E^ 



1 1 
1 1 



and Fo = 



Ej I n 



iT 



. From the conditions of Theorem 4.1 we find the piece- 



wise quadratic Lyapunov function V(x) = x T PiX with 



Pi = Ps 



5 0 
0 1 5 



P2=P4 



1 0 
0 5 



As seen in Figure 4. 7 (right), the level surfaces of the computed Lyapunov 
function are neatly tailored to the system trajectories. □ 



Analysis of a Min-Max Selector System 

The examples analyzed so far have been small examples in two dimensions, 
constructed to illustrate the shortcomings of quadratic and merits of piece- 
wise quadratic Lyapunov functions. Our next example is motivated by in- 
dustrial applications, has higher state dimension (seven continuous states), 
and a nonlinearity that is not easily dealt with using alternative techniques. 
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Figure 4.8 Control system with min/ max selectors, from [8] . 



The example is the min-max selector control system shown in Figure 4.8. 
This scheme is common in situations where several process variables have 
to be taken into account using a single control signal. In Figure 4.8, y is the 
primary variable and z is a process variable that must remain within given 
ranges. The controller C is designed to control the primary variable while 
the controllers C max and C min are designed to keep the critical variable z 
within certain bounds. Designed correctly, the min-max selector chooses the 
controller that is most appropriate for the situation and allows good control 
of the primary variable while respecting the constraints. 

Consider a system characterized by 



GT(s) 

G2(S) 



40 

0.05s 3 + 2s 2 + 22s + 40’ 
5 

s 2 + 7s + 5 



To control the primary variable, we design a lead-lag controller 



s 2 T 3s T 3 
= 0.02s 2 + s + 0.01 7 



while Cjnin and C max are proportional controllers 



Ul = Ki(z m iri - z). 
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The plots in Figure 4.9 show a simulation of the system without constraint 
controllers. The tracking of the primary variable is quite good, but the criti- 
cal variable 2 exceeds it constraint limits (shown in dashed lines). The plots 
in Figure 4.10 show a simulation of the min-max selector strategy. The con- 
straints are now respected while the tracking of the primary variable re- 
mains satisfactory. 




0 10 20 30 40 



Figure 4.9 Simulation of the control system in Figure 4.8 without constraint handling. 

The tracking of the primary variable is quite good (top) , but the critical variable ex- 
ceeds it limits (bottom) . 

We will apply Theorem 4.1 to stability analysis of the system for a con- 
stant set-point y sp and constant constraint limits z m in , z ma x- Different values 
of y sp , Zmin and z m ax result in different equilibrium points. For sake of sim- 
plicity, we will let y sp = z m ax = Zmin = 0, but the technique would apply 
similarly to any choice of constant inputs. 

For analysis purposes, it is convenient to re-write the system equations 
as a linear system interconnected with the static nonlinearity 

u = min (uh, max (w n , u {)) . 

The selector nonlinearity has three input signals u n , ui and one output, u. 
Similarly to the selector system used as motivating example in Section 4.4, 
we can reduce the number of inputs by one using a simple loop transfor- 
mation. This results in the system shown in Figure 4.11. The transformed 
system has two outputs, vm and v n i, and the selector nonlinearity is now 
reduced to the two-dimensional mapping tp(vhi,Vni) shown in Figure 4.12. 
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0 10 20 30 40 



Figure 4.10 Simulation of the min-max selector system. The constraint limits are re- 
spected (bottom) , while the tracking of the primary variable is still satisfactory (top) . 
The middle plot shows how the constraint controllers override the primary control 
signal (dashed) resulting in a control (full) that respects the constraints. 




Figure 4.11 Selector control system rewritten as linear system interconnected with a 
static multi-variable nonlinearity. 



The nonlinearity p is piecewise linear and has the explicit expression 

f 0 if (vhi > 0) A (v n i < 0) 

p = \ v n i if ( v n i > 0) A ( Vhi > v n i ) 

[ vm otherwise 

Since the region where <p(vhi,v n i) = Vhi is not convex, we have to introduce 
an additional region, see Figure 4. 12 (right). 

While this nonlinearity fits directly in the piecewise linear framework, it 
is not easily dealt with using other techniques. It is easy to verify that the 
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Figure 4.12 Static nonlinearity (left). The corresponding non-convex state partition 
(right) is rendered convex by splitting one cell in two (the dashed line in rightmost 
figure) . 



nonlinearity has gain less than one, which motivates an attempt to apply 
the small gain theorem. However, the £ 2 mduced gain of the linear system 
is 15.8, and the small gain can not verify stability An approach based on 
linear differential inclusions (Corollary 4.1) also fails. 

In contrast, a numerical stability analysis using Theorem 4.1 verifies sys- 
tem stability. In this case, the optimization returns a Lyapunov function 
which is globally quadratic. Since Corollary 4.1 fails to establish stability, 
this example shows the importance of using partition information in the 
analysis. 



4.7 Analysis of Piecewise Linear Differential Inclusions 

The LMI computations of Theorem 4.1 are readily extended to systems de- 
scribed by piecewise linear differential inclusions (pwLDIs) , 

x(t) € co {. Akxit ) + dk} x(t) € Xi. 

keK(i) 

In order to guarantee that the Lyapunov function is decreasing along all pos- 
sible solutions of the pwLDI, we must require that the Lyapunov function is 
decreasing with respect to every extreme dynamics x = AkX + a* that de- 
fines the inclusion in each cell. This leads to multiple decreasing conditions 
in each region (one for each extreme dynamics) and the following result. 

Theorem 4.2— PwQ Stability of PwLDIs 

Consider symmetric matrices T, U \k and such that U \k and have 
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Figure 4.13 The left figure shows simulations (full) and Lyapunov function level sur- 
faces (dashed) obtained in Example 4.8. The right figure shows the computed Lya- 
punov function. 



nonnegative entries, while 

Pi = Fj TFi , i £ Jo 

Pi = PJ TFi > i £ I\ 



satisfy 



0 > Aj^Pi + PiAk + Ej UikEi 
0 <Pi -EjW ik Ei 

0 > Ak Pi + PiAk + Ej UikEi 

0 <Pi -EjW ik Ei 



i £ Iq , & G iT(i) 



i £ -Zi, f G K(%) 



Then every solution x(t) G of the inclusion (2.12) with u = 0 for t > 0 

tends to zero exponentially □ 

Proof Follows similarly to Theorem 4.1 and Proposition 4.2. 



Theorem 4.2 enables the use of the piecewise quadratic machinery for 
rigorous stability analysis of smooth nonlinear systems. The following ex- 
ample demonstrates the ideas. 

Example 4.8— Embedding Smooth Systems in PwLDIs 
S imulations indicate that the following nonlinear system is stable 
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x\ = — 2xi + 2 x 2 + sat(xiX 2 )xi 
x 2 = ~2xi - sat(xix 2 )(xi + 4 x 2 ). 





4.8 Analysis of Systems with Attractive Sliding Modes 

We would like to verify global exponential stability of the origin by com- 
puting a piecewise quadratic Lyapunov function for the system. A simple 
technique for rigorous analysis of the system is to explore the bounds 



Pmin L Sat (x 1^2) Li Pm ax 



and re-write the model as the differential inclusion 



x = 



-2 

-2 



2 

0 



x+p(t) 



1 

-1 




(4.25) 



with Pmin < p(t) < Pmax • By using information about the nonlinearity 
p(t) = sat(xiX 2 ) we can obtain pwLDIs of different accuracy. First, notice 
that analysis using a global model based on the bound — 1 < p(t) < 1 is 
futile since p(t) = -1 gives an unstable extreme dynamics. Taking the step 
from linear analysis to piecewise linear analysis, we can obtain a refined 
model by exploring the fact that 

0 < sat(xiX2) < 1 

in the first and third quadrant, and 

—1 < sat (x 1 x 2 ) < 0 

in the second and fourth quadrant. This observation motivates a model with 
four regions, each region covering one quadrant. The dynamics in each re- 
gion is given by a linear differential inclusion on the form (2.12). To verify 
stability, we apply Theorem 4.2 and find the Lyapunov function with the 
level curves indicated in Figure 4.8. This proves global exponential stability. 
Note that the level surfaces are non-convex sets and that the system is not 
easily dealt with using absolute stability results due to the multi-variable 
nature of the nonlinearity sat(xiX 2 ). □ 



4.8 Analysis of Systems with Attractive Sliding Modes 

The analysis procedures developed so far are only concerned with systems 
that do not have any attractive sliding modes. As demonstrated by the fol- 
lowing example, however, attractive sliding modes can play a critical role 
for system stability and must be accounted for in the analysis. 
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Example 4.9— Sliding Modes and Stability 
Consider again the piecewise linear system from Example 2.4: 



A\x 



x = < 



A 2 x 




-2 

-4 




x\ > 0 



x\ < 0 



(4.26) 



Since both system matrices are stable, and the system admits a solution to 
the LMIs in Theorem 4.1, it is easy to be misled to believe that this would 
imply exponential stability of (4.26). However, we know from Example 2.5 
that the system does not have any non trivial trajectories, so Theorem 4.1 
does not apply. In fact, we also know that the dynamics on the sliding mode 
is unstable, x 2 = x 2 , so once Filippov solutions reach the attractive sliding 
mode, they tend to infinity along the positive x 2 -axis, see Figure 4.14. □ 




Figure 4.14 The behavior on attractive sliding modes is critical for system stability. 
The Filippov solution reaches the sliding mode and tends to infinity. 



For systems with attractive sliding modes, it is natural to look for conditions 
that ensure stability of Filippov solutions. To this end, let 



5 = {x | Gx h 0 A Hx = 0} (4.27) 

be a set on which there is an attractive sliding mode (this might be a regular 
sliding set on the boundary between two cells, a higher-order sliding set, or 
a sliding set on the intersection of multiple boundaries) and define 

K s = {k\x k ns^<b}. 
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4. 8 Analysis of Systems with Attractive Sliding Modes 
Thus, any Filippov solution that remains in S satisfies 

x(t) e co {A k x{t) + a k } 

k^K s 

Let V (x, i) = x T PiX be the expression used for defining the Lyapunov func- 
tion candidate in cell Xi. If, for some i e K s , 

dV g*’^ {A k x + a k }< 0 Vx e <S\{0}, \/k€K s (4.28) 

then V ( x ) is decreasing along Filippov solutions on S . This condition can be 
reformulated as the LMI 

AjPi + PiA k + G T s W ik G s + Hj N ik + Nj k H s < 0 \/k e K. (4.29) 

where W ik >z 0 and N ik are matrix variables of appropriate dimensions. 
Thus, the piecewise quadratic analysis can be extended to systems with at- 
tractive sliding modes by first identifying sliding sets on the form (4.27) and 
then augmenting the conditions of Theorem 4.1 with constraint of the type 
(4.29). Due to continuity of V(%), we only need to enforce the above condi- 
tions for one i e K s (if V(x,i) is decreasing then so is V (x, j ) for all j € K s ) . 
The following example illustrates the approach. 

Example 4.10— Exponential Stability of Sliding Mode System 
Consider the following slight modification of Example 4.9 



A\x 


X2 > 0 


A 2 X 


X2 < 0 



Here, and A 2 are the system matrices defined in Example 4.9. Also this 
system has an attractive sliding mode as illustrated by the Filippov solution 
shown in Figure 4.15 (left). To verify exponential stability of all Filippov 
solutions, we apply the above technique and find the LMI conditions 



Aj 1 P 2 + P^A\ + 



0 



A^Pi -\~ P\A\ < 0 
Al 2 P 2 + P 2 A 2 < 0 



l] T Ny N 



0 



1 < 0 



Pi >0 
P2 > 0 



where Pi = FjTFi for i = 1,2. The conditions have a feasible solution 
which guarantees that all Filippov solutions converge exponentially to zero. 
Time responses of a Filippov solution and the corresponding value of the 
computed Lyapunov function are shown in Figure 4.15. □ 
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Figure 4.15 The extended stability conditions verify stability also in the presence of 
attractive sliding modes. The computed Lyapunov function decreases also on the at- 
tractive sliding mode (right) . 



Attractive sliding modes give a significant increase in the analysis compu- 
tations since we need to analyze a differential inclusion on every sliding set, 
and sliding sets can potentially appear on every cell boundary and inter- 
section of boundaries. Hence, if sliding modes can be ruled out a priori one 
should use theorems that only considers trajectories. For systems with slid- 
ing modes it is advised to first rule out attractive sliding modes on as many 
boundaries as possible before applying the methodology described above. 



4.9 Improving Computational Efficiency 

We have seen how the piecewise quadratic analysis is much more powerful 
than the classical approach based on quadratic Lyapunov functions and how 
it allows us to analyze many systems where other methods fail or are hard to 
apply. However, the piecewise quadratic approach is computationally more 
demanding and a straightforward implementation of the LMI conditions in 
Theorem 4.1 may result in very large optimization problems (this is espe- 
cially true when the state space partitioning is performed in many dimen- 
sions). It is therefore of interest to develop methods that decrease the com- 
putational burden without introducing excessive conservatism. This section 
will describe two such approaches. 

Stability Analysis in Two Steps 

The LMI conditions in Theorem 4.1 incorporate the positivity condition in 
the Lyapunov function search. Looking back at Lemma 4.1, however, we see 
that there is little reason for doing so. The result suggests that any function 
which is decreasing along system trajectories contain all information about 
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system stability. If the function is positive definite, stability follows analo- 
gously to Theorem 4.1. If we find some point where the computed function 
is non-positive, on the other hand, then no trajectory passing through this 
point can approach the origin as t — > oo. We give the following result. 

Proposition 4.4— Stability Analysis in Two Steps 
Consider a symmetric matrix T and symmetric matrices Ui y 0 with non- 
negative entries, while Pi = FfTFi and Pi = FjTFi satisfy 

0 >AjPi + PiAi + Ej 1 UiEi i G Jo 

0 >Aj Pi + PiAi + Ef UiEi i € /i 

Let x(t) G U i e iXi be a trajectory of (2.3) with u = 0 for t > 0, and define 

V (x) = x T PiX for x G Xi, i G /. 

If there exists a > 0 such that V{x) > a||x||2 , then every trajectory x(t) tends 
to zero exponentially. If V(xo) < 0 for some xo G X\{0}, then no trajectory 
x(t) with x(0) = xo can tend to zero as t — ► oo. □ 

Proof: Follows from Lemma 4.1 along the lines of the proof of Theorem 4.1. 

Proposition 4.4 implies that the large LMI problem in Theorem 4.1 can 
be split into several smaller problems. By disregarding the positivity con- 
straints in the Lyapunov function search we reduce the number of con- 
straints by roughly 50% and eliminate the need for the variables W it This 
problem can be solved in a fraction of the time needed to solve the original 
problem. Moreover, if the LMI conditions in Proposition 4.4 do not admit a 
solution, then neither do the analysis conditions in Theorem 4.1. 

Once a Lyapunov function candidate is found we can use Lemma 4.3 to 
check if this function is positive on the partition. Verifying positivity then 
amounts to solving a single LMI problem for each region. Since the Lya- 
punov function is now fixed, each such problem has only one constraint in 
one free matrix variable and can be solved very efficiently. We will illustrate 
the savings obtained by Proposition 4.4 in the end of this section. 

Quadratic Cell Boundings - Computational Savings at a Price 

It is often the number of free parameters in the 5-procedure relaxation (the 
entries of the matrices Ui and Wf) that adds the most parameters to the Lya- 
punov function search. A second approach for reducing the computational 
burden is therefore to limit the number of parameters in the 5-procedure. 
Returning to Lemma 4.3 we see that for given Pi a solution to the inequality 

Pi - EfUiEi > 0 
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implies that V ( x ) is positive for all x in the quadratic set 
Si = {x | x T Ej UiEix > 0}. 

Since Xi C £ it we may view Si as a quadratic cell bounding, see Figure 4.16. 
From this perspective, the free parameters in Ui are used to form a quadratic 
set that separates Xi from the set V~ = {x \V(x) < 0}. One way to reduce 
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Figure 4.16 Several quadratic boundings £i(Ei) (dark) of the cell Xi (light) can be 
derived by varying the entries in the matrix U i . 

the number of search variables is to fix a quadratic cell bounding before 
the Lyapunov function search. Procedures for computing optimal ellipsoidal 
approximations of polyhedra are described in Appendix A. To pursue this 
direction further, we define quadratic cell boundings as follows. 

Definition 4.3— Quadratic Cett Bounding 
A matrix Si = Sj is a quadratic cell bounding if 

x T Six > 0 for x(t) G Xi. (4.30) 

Furthermore, we say that {*%} have the zero interpolation property if 



h nx l 

Olxn 0 



for i G /q. (4.31) 



□ 

Now, rather than using Lemma 4.3 in the conditional analysis we use the 
following result. 



Lemma 4.4 

Consider a piecewise quadratic function 



V(x) 



' X T PiX 



C T Pi X 
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x G Xi i G Io 
x G Xi i G I\ 



(4.32) 
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with Pi = Pj, Pi = PT and let Si be a quadratic cell bounding satisfying 
(4.30) and (4.31). If there exists positive scalars Wi > 0 such that 



Pi - wA > 0 (4.33) 

Pi - wA > 0 (4.34) 

then there exists a > 0 such that V (x) > a\\x\\\ lov d\\ x e X . □ 

Proof: Follows similarly to Lemma 4.3. 

The following variant of Theorem 4.1 now follows directly. 



Proposition 4.5— PwQ Stability with Quadratic Relaxation 
Consider a symmetric matrix T and nonnegatvive scalars Ui and Wi such 
that Pi = Fj'TFi for i e I 0 and Pi = F^TFi for i e h satisfy 



0 > Aj Pi + PiAi + UiSi 
0 <C Pi WiSi 

0 > Af Pi + PiAi + Ui§i 
0 <C Pi WiSi 



i e Io 



i e h 



Then every trajectory x(t) G U ie iXi satisfying (2.3) with u = 0 for t > 0 
tends to zero exponentially. □ 

Proof: Follows similarly to Theorem 4.1. 

This approach allows for large savings search variables, but is more con- 
servative than the original formulation in Theorem 4.1. To be more precise, 
assume that Ei e RP x ( n + 1 ), Then the polyhedral relaxations EfUiEi use 
p(p — l)/2 free parameters, while the quadratic formulation (4.30) uses one 
single parameter. The conservatism comes from the requirement that the a 
quadratic cell bounding has to be fixed before the optimization. 

A natural candidate for quadratic approximation of polyhedral cells is to 
compute the ellipsoid with minimum volume that contains the cell. How- 
ever, minimal volume has little to do with the role of the relaxation term 
in the LMI conditions, and in Chapter A we will be able to prove that for 
some important classes of partitions the use of minimal volume ellipsoids 
and Proposition 4.5 is always more conservative than the original formula- 
tion in Theorem 4.1. Such a proof requires some further developments, but 
we can already now demonstrate the arguments on a simple example. 
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A Comparative Example 

To give a flavor of the benefits and limitations of the different formulations 
of piecewise quadratic stability given in Theorem 4.1, Proposition 4.4 and 
Proposition 4.5, we will apply them to the system shown in Figure 4.17. The 
system dynamics is given by 



x = Ax + 6i/i(xi) + 62/2(2:2) 

where A e R 2x2 , 61,62 £ R 2xl and fi(xi) = arctan (x*). We will present 





Figure 4.17 The system used as comparative example (left). The nonlinearity fi(xi) 
is shown in full lines in the right figure. The dash-dotted line illustrates a piecewise 
linear approximation and the dashed lines show piecewise linear sector bounds. 

results for both piecewise linear approximations and piecewise linear sector 
bounds on the nonlinearities, see Figure 4.17 (right). In both cases, the piece- 
wise linear descriptions induce a partition of the domain [-5, 5] x [-5, 5] into 
nine regions. 

First, we let 




and use the piecewise linear approximation of fi(xi). In this case, all ap- 
proaches succeed in verifying stability. The computational results are shown 
in Table 4.1. The computations were performed on a SUN Ultra 10 com- 
puter using the LMI software [62]. In the table the acronym P refers to the 
use of polytopic cell boundings in the analysis while Q indicates the use of 
quadratic cell boundings. The number 1 means that the analysis was per- 
formed in a single step (enforcing both positivity and decreasing conditions 
simultaneously) while 2 means that the analysis was performed in two steps 
(enforcing the decreasing condition during the Lyapunov function search 
and subsequently verifying positivity). 
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Table 4.1 First set-up: all approaches verify stability. Large savings in computations 
are obtained from the alternative formulations (P-2, Q-l, Q-2). 



Approach 


Time (s) 


^Variables 


# Constraints 


p-i 


1.04 


117 


114 


P-2 


0.41 


69 


57 


Q-l 


0.23 


37 


34 


Q-2 


0.11 


29 


17 



As seen in Table 4.1, Proposition 4.4 (P-2) results in a large reduction 
in computation time compared with the computations required by Theo- 
rem 4.1 (P-1). The computational savings are even greater when using quad- 
ratic cell boundings as in Proposition 4.5 (Q-l). In this case, the quadratic cell 
boundings are taken as the minimal volume ellipsoids that cover each region 
(see Proposition A. 10). By combining the two-step analysis with quadratic 
cell boundings (Q-2) the computational time is reduced to around 10% of 
what was required by the original formulation. 

Using the same matrices A : b\ and & 2 > we now consider the case when 
the nonlinearities are described by piecewise linear sector bounds. This ap- 
proach allows us to verify stability of the smooth nonlinear system in a rig- 
orous way but it also increases the computational cost. In each region the 
system is now described by a differential inclusion with four extreme dy- 
namics. As the main burden in analysis of such systems is verification of the 
multiple decreasing conditions, the savings of the two-step analysis proce- 
dure is somewhat reduced, see Table 4.2. 



Table 4.2 Second set-up: the use of piecewise linear sector bounds decreases the ben- 
efits of the two-set analysis procedure, but good savings are still obtained. 



Approach 


Time (s) 


#Variables 


# Constraints 


p-i 


3.79 


261 


285 


P-2 


2.17 


213 


228 


Q-l 


0.80 


61 


85 


Q-2 


0.45 


53 


58 



The drawback with quadratic cell boundings is that they allow very little 
freedom in the 5-procedure-relaxation. This introduces some conservatism 
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as can be seen by letting 




and using piecewise linear sector bounds on the nonlinearities. The compu- 
tational results are shown in Figure 4.3. Stability can no longer be verified 
using quadratic cell boundings, while the computational savings in the use 
of Proposition 4.4 remain the same. 



Table 4.3 Final set-up: quadratic cell boundings fail to verify stability. 



Approach 


Time (s) 


#Variables 


# Constraints 


p-i 


4.15 


261 


285 


P-2 


2.62 


213 


228 


Q-l 


fails 


- 


- 


Q-2 


fails 


- 


- 



To understand the computational complexity of the different approach- 
es, it is useful to see how different factors contribute to the total number 
of parameters. In this example, we have constructed the constraint matrices 
using the procedure given in Section A.4. This procedure gives Fi e R 6x3 
and Ei G R 4x3 . This implies that T e R 6x6 and that the Lyapunov func- 
tion candidate F^TFi has 21 free parameters. Each of the matrices Ui and 
Wi used in the polytopic 5-procedure has 6 free parameters while the ellip- 
soidal 5-procedure uses 1 parameter. As the origin lies in the interior of one 
cell, 5-procedure relaxation is only used in 8 regions. Applied to the first set- 
up (Table 4.1), the approach P-1 requires 21 + 8 * (1 + 1) * 6 = 117 parameters 
while P-2 uses 21 + 8 • 1 • 6 = 69 parameters. For the piecewise linear sector 
bounds (Table 4.2), P-1 uses 21 + 8 • (1 + 4) • 6 = 261 parameters while P-2 
uses 21 + 8-4-6 = 213 parameters. In this case, Q-2 uses only 21 + 8-4 - 1 = 53 
parameters. 



4.10 Piecewise Linear Lyapunov Functions 

There are several reasons to look for alternatives to the piecewise quadratic 
analysis. First, the semidefinite programming problem in Theorem 4.1 has 
many free search variables and analysis of large problems using these con- 
ditions may be time consuming using today’s LMI software. Secondly, the 
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5-procedure, which was used to exploit the restriction x e Xi in the compu- 
tations, can be a restrictive way to check positivity of a piecewise quadratic 
function on a polyhedral partition. In other words, there are systems that 
admit a piecewise quadratic Lyapunov function but where this function can 
not be found by the formulation of Theorem 4.1. 

As an alternative to piecewise quadratics, we will consider Lyapunov 
function candidates that are continuous and piecewise linear 



V(x) 



pfx x e Xi , i G io 

pfx = pj x -\- qi x e Xi, i £ I 1 



(4.35) 



Such Lyapunov functions can be computed via linear programming. Com- 
pared to LMI software, linear programming solvers have reached a high 
level of maturity and efficient software exists for large-scale problems that 
allows systems with thousands of cells to be analyzed in a matter of seconds. 

Similar to the piecewise quadratic case we will use a compact param- 
eterization of such functions that separates the free parameters from the 
constraints imposed by the continuity requirement. The parameterization 
is established in the following lemma. 



Lemma 4.5— PwL Parameterization 

Let { Xi }i € i be a polyhedral partition, and let Fi € R^ x ( n+1 ) be continuity 
matrices that satisfy (4.14). Then, for each t e R p , the scalar function 

V (x) = t T Fix := piX for x e Xi. 

is continuous and piecewise linear. Moreover, if {Fi} have the zero interpo- 
lation property, there exists a and (3 such that 

Moo < V(x) < (3 Moc- 



□ 

Proof: Follows similarly to Lemma 4.2 and the absence of offset terms in V (x) 
in a neighborhood around the origin. 

Another attractive feature of piecewise linear Lyapunov functions is that 
the conditional analysis can sometimes be performed without loss. This fact 
is established in the following Lemma, similar in nature to Farkas’ lemma 
[181,221]. 



Lemma 4.6 

The following statements are equivalent. 
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1 . p T x > 0 for all x such that Ex y 0, Ex ^ 0 

2. There exists a vector u >- 0 such thatp - E T u = 0. 

□ 



Proof: See Section B.2. 

Contrary to the standard formulation of Farkas’ lemma which uses non- 
strict inequalities, Lemma 4.6 is formulated using strict inequalities. How- 
ever, the result only considers linear forms and does not treat affine terms. 
The following stability theorem now follows. 



Theorem 4.3— Piecewise Linear Stability 

Let {Xi} ie i C R n be a polyhedral partition with continuity matrices F i} sat- 
isfying (4.14) and (4.15), and cell boundings E i} satisfying (4.19) and (4.20). 
Assume furthermore that Erx ^ 0 for every xgI* with x ^ 0. If there exists 
a vector t and non-negative vectors Ui y 0 and wi>- 0 while 



Pi = F?t, 


i € Io 




Pi = Fft, 


i e h 




j 0 = pf A{ + iiiEi 
\o =pj - WiEi 


i e I 0 


(4.36) 


j 0 = pf Ai + UiEi 
[o = pj - WiEi 


i e h 


(4.37) 



then every trajectory x(t) e U ieI Xi satisfying (2.3) with u = 0 for t > 0 tends 
to zero exponentially. □ 

Proof: Follows similarly to Theorem 4.1. 

The search for free variables t, ui and in Theorem 4.3 is a linear pro- 
gramming problem. If the system does not have any attractive sliding modes 
a solution to this problem guarantees that 



V(x)=pfx for xeXi.iel 

is a Lyapunov function for the system. Although the analysis conditions still 
use relaxation terms ui and w it the number of entries in ui and wi has been 
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reduced in comparison to the number of entries in the matrices Ui and Wi 
used in the piecewise quadratic analysis. Finally, note that the additional 
constraint 



EiX ^ 0 for x e Xi with x ^ 0 

does not impose any further restriction. For i e Io, the assumption is vio- 
lated if {x | EiX y 0} is some linear halfspace, but p T x can not be strictly 
positive for all x in a closed linear halfspace. For i g Ii, the situation can 
always be avoided by adding the additional constraint [0i xn l]x > 0 to 
the cell boundings. 

Piecewise Linear Lyapunov Functions on Bounded Partitions 

Theorem 4.3 can be used for analysis of systems with both bounded and 
unbounded partitions. However, if all cells are bounded we can reduce the 
computations even further. More specifically, let the cells be given in vertex 
representation 



Xi = CO {is k } 

kev(i ) 

where V (i) are the set of indices for the vertices v k of cell X it Then, an affine 
function is positive on Xi if and only if it is positive on the vertices of X 
This allows us to formulate the following result. 

Theorem 4.4 

Let {Xi} ie i be a partition of a bounded subset of R n into convex polytopes 
with vertices i/ k , and let Fi be the associated continuity matrices satisfying 
(4.14) and (4.15). If there exists a vector t such that 

Pi = Fjt for % G Io 

pi = Fjt for i e I 



satisfy 



> pf An/, k 


i G /o, 


i/ k G Xi 


<pjvk 


i £ Iq-> 


i/ k G Xi 


> pJAiPk 


i £ lu 


v k G Xi 


<pjv k 


i G /i, 


is k G Xi 



(4.38) 



(4.39) 
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for each v ^ 0, then every trajectory x{t) e U ie iXi satisfying (2.3) with 
u = 0 for t > 0 tends to zero exponentially □ 

Proof: Follows similarly to Theorem 4.3 but where decreasing and positivity 
conditions are checked according to the discussion above. 

Note that all the relaxation terms have vanished, and that the vector in- 
equalities of Theorem 4.3 have been reduced to scalar inequalities. 

It is possible to arrive at even simpler stability conditions if one considers 
partitions where each cell Xi C R n has n + 1 vertices. Such polytopes are 
called simplices, and are described in more detail in Chapter A. For simplex 
partitions, the Lyapunov function is completely determined by its values at 
the cell vertices, and the positivity condition can then be replaced by the 
requirement that all entries of the vector t should be positive, t >- 0. 

Polytopic Lyapunov Functions 

Although Theorem 4.4 requires the analysis domain to be bounded (the cells 
are polytopes) it can still be used to assess global exponential stability in 
some cases. More precisely, if A = 0, then any piecewise linear Lyapunov 
function valid in some open neighborhood of the origin is also valid glob- 
ally. Lyapunov functions derived in this way are often called polytopic Lya- 
punov functions, as the level sets of such a Lyapunov function are polytopes 
(see [23, 147] for further details). The following example illustrates the ideas. 

Example 4.1 1— Selector System Cont’d 

To illustrate the use of piecewise linear Lyapunov functions, we return to 
the simple selector system used in Section 4.4. As discussed in conjunction 
with Theorem 4.3, the piecewise linear Lyapunov functions cannot be used 
on the initial partition since the natural cells are both open linear halfspaces. 
Using the refined partition shown in Figure 4.18 (left) , however, Theorem 4.4 
returns the Lyapunov function shown in Figure 4.18 (right). Hence, global 
exponential stability follows from the arguments above. Note that the Lya- 
punov function is poorly conditioned and that the level surfaces are heav- 
ily unbalanced. By refining the partitioning further one arrives at Lyapunov 
functions that closely resemble the Lyapunov function used in the piecewise 
quadratic analysis (not shown). □ 



Extensions 

The basic stability computations can be extended in several useful ways. 
One example is computation of decay rate, r, which can be estimated from 
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Figure 4.18 Refined partition and level surfaces of computed Lyapunov function 
(dashed) to the left. The computed Lyapunov function is shown to the right. 



the modified Lyapunov inequality 

V (x) + rV ( x ) <0 Vx 0. 

Given a fixed value of r, the above condition can be verified using a slight 
modification of the previous theorems (where Ai has been replaced by Ai + 
rl in the decreasing conditions). The optimal value of r can then be found by 
bisection. Another possibility is to extend the stability analysis to piecewise 
linear inclusions, 

x <G co { A k x-\-a k } x e Xi 

keK(i) 

In this case, one simply needs to verify multiple decreasing conditions in 
each region. Similar techniques can be used to assess stability of systems 
with attractive sliding modes. 



4.11 A Unifying View 

There is a close relationship between the parameterizations of the piecewise 
linear and the piecewise quadratic Lyapunov functions. In this section, we 
will elaborate this relationship further and establish a unifying framework 
for computation of globally quadratic, piecewise quadratic, polyhedral and 
piecewise linear Lyapunov functions. 

One may view the matrix format for continuous piecewise quadratic 
Lyapunov functions introduced in Lemma 4.2, 

V(x) = x T FjTFiX x € 
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as a quadratic form in the coordinates z obtained by the continuous piece- 
wise linear mapping z = F{X. If one rather considers linear forms in z, 

V{x) = t T z = t r FiX := pjx x E Xi 

one obtains the parameterization of continuous piecewise linear functions 
suggested in Lemma 4.5* In fact, the piecewise linear Lyapunov functions 
can be seen as a direct restriction of the piecewise quadratics as follows. Let 
!*\ be constraint matrices satisfying (4*14) and define 



T si 




0 


t 


0 1 


: 2 


t T 


0 



Note that P* also satisfies the continuity condition (4*14), and that the piece- 
wise quadratic function V(x) of Proposition 1 now evaluates to 

x T P l x = ^FfTFiX = t T FiX = p[x. 

The close relationship between the parameterizations of piecewise lin- 
ear and piecewise quadratic Lyapunov functions allows us to establish a 
unifying view of several approaches for numerical Lyapunov function con- 
struction, see Figure 4.19. 



V(x) = aFFTTFiX 




Figure 4.19 A unifying view of four classes of Lyapunov functions for piecewise lin- 
ear system that can be computed via convex optimization. 

Most versatile are the piecewise quadratic Lyapunov functions [100, 157, 72] 
V[x) = ;c T FTTFrx x € Xi 
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As shown in Theorem 4.1, piecewise quadratic Lyapunov functions can be 
computed via convex optimization in terms of LMIs. The fact that a certain 
set of differential equations only describes the dynamics in a restricted re- 
gion of the state space can be exploited using the 5-procedure, which is a 
sufficient condition that appears to work very well in practice. 

The quadratic Lyapunov functions [45, 27] are special instances of the piece- 
wise quadratics obtained by letting Fi = [/ nxn On]- Quadratic Lyapunov 
functions can be computed via LMI optimization, and conditional analysis 
can be done using the 5-procedure. 

Also the piecewise linear Lyapunov functions [118, 96, 107] 

V(x)=t T FiX xGli 

can be seen as a special case of the piecewise quadratics. They can be com- 
puted via linear programming as shown in Theorem 4.3 and Theorem 4.4. 
The conditional analysis can in some cases be formulated without loss (as 
established in Lemma 4.6). 

Polytopic Lyapunov functions [142, 147, 23] are a special case of the piece- 
wise linear Lyapunov functions. The polytopic Lyapunov functions can be 
obtained from the piecewise linear by considering partitions that consist of 
convex cones with base in the origin, i.e., polyhedral partitions for which 
Ji = 0. The computations can be done using linear programming and the 
conditional analysis is performed without loss. 

Choosing Lyapunov Function Class 

The choice of Lyapunov function class involves many trade-offs between 
issues such as analytic flexibility, complexity of description, and computa- 
tional requirements. For example, the quadratic Lyapunov functions have 
a compact description and allow for efficient computations but has limited 
analytic flexibility. Piecewise linear and piecewise quadratic Lyapunov func- 
tions have a high degree of flexibility but are more demanding to compute 
and require more parameters for their representation. 

For a given partition, the linear matrix inequalities in Section 4.6 are 
much more demanding to solve than the linear programming problems in 
Section 4.10. However, as we have seen, it is often necessary to refine an 
initial partition in order to prove stability using a piecewise linear Lya- 
punov function. Partition refinements increase both the computational re- 
quirements and the number of parameters needed for representing the Lya- 
punov function. A particular weakness of piecewise linear Lyapunov func- 
tions is that a very fine partition is needed when the system dynamics is 
oscillatory. The following example illustrates this issue. 
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Example 4.12— Lyapunov Functions and Oscillatory Dynamics 
Consider the linear oscillator 



—a f3 
(3 —a 



(4.40) 



The system matrix has eigenvalues A = -a ± i(3 and the system is expo- 
nentially stable (hence, admits a quadratic Lyapunov function) for all a > 0. 
The larger the value of /?, the more oscillatory the dynamics and the more 
circular the trajectories in the phase plane. Circular trajectories require a Lya- 
punov function with circular level sets, and if we want to approximate these 
by polytopes (which is what we do when we use polytopic Lyapunov func- 
tions) the number of polygon sides might need to be very large [166]. In fact, 
Polanski [161] proved that a poly topic Lyapunov function with a partition 
defined by the 2 m vertices 



± 




( cos(7r/ra) 
sin(7r/m) 



( cos(7r(m — 1 )/m) 
sin(7r(m — 1 )/m) 



is a Lyapunov function for the system (4.40) if 



7 r 

'TY) 

arctan (2 a//?/ (1 — a 2 / (3 2 )) 

The number of vertices grows large when a/ f3 grows small, see Figure 4.20 
and tends to infinity as a/ j3 tends to zero. □ 




Figure 4.20 Polanski’s construction for a//3 = 1 (left) requires 6 vertices (marked 
with 'x'), a//3 = 0.1 requires 32, while a//3 = 0.01 (right) requires 316 vertices. All 
three systems admit a quadratic Lyapunov function that can be described by three 
parameters, and computed by solving a single Lyapunov equation. 



The main reason for moving beyond quadratic Lyapunov functions is to in- 
crease flexibility of the Lyapunov function candidate. This flexibility is not 
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needed for analysis of linear systems but it can be critical in the analysis of 
nonlinear and uncertain systems. While the level sets of quadratic Lyapunov 
functions are always ellipsoidal, poly topic Lyapunov functions allow level 
sets of arbitrary convex polyhedral shape. The piecewise linear and piece- 
wise quadratic Lyapunov functions are even more flexible since they can be 
both non-convex and non-homogeneous. The following example, also used 
in [150, 23], demonstrates some of these issues. 

Example 4.13— Analytic Flexibility 
Consider the following linear uncertain system 



x(t) = A{S(t)}x(t) 



0 

— 1 + 6(t) 



1 

-1 



x{t). 



where |<5(f)| < d is an uncertain time-varying parameter. For d > a/3/2 there 
is no quadratic Lyapunov function that can show stability for all admissi- 
ble parameter variations. In [23] , Blanchini reported a polyhedral Lyapunov 
function with 30 vertices that proves stability for d < 0.98. Using a piecewise 
quadratic Lyapunov function with four regions (being the four quadrants in 
R 2 ) and Theorem 4.2 we can not only decrease the number of parameters 
needed to represent the Lyapunov function but also improve the bound to 
d = 1 — e with e = IE — 15. □ 

We believe that the piecewise quadratic Lyapunov functions strike a nice 
tradeoff between analytic flexibility and complexity of description: they are 
very versatile (as they encompass all other Lyapunov function classes dis- 
cussed in this book) and they allow analysis of systems with oscillatory dy- 
namics without the need for excessive partition refinements. However, they 
are computationally more demanding and the S-procedure introduces many 
parameters in the optimization problem. If a very fine partition is needed in 
order to approximate an underlying nonlinear system by a piecewise linear 
system and the linearized dynamics around the equilibrium point is not too 
oscillatory, piecewise linear Lyapunov functions might be the best choice 
(see also the Notes and References section at the end of this Chapter). 

Another issue appears when Lyapunov-like functions are used in system 
analysis and optimal control problems. Different problem formulations then 
call for different classes of loss functions. While energy-related problems, 
such as computation of the induced £2 -gain, are conveniently dealt with us- 
ing quadratric or piecewise quadratic Lyapunov functions, piecewise linear 
and polytopic Lyapunov functions have been useful in analysis of systems 
with absolute constraints (c.f., Section 3.2). 
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4.12 Comments and References 



Numerical Lyapunov Function Construction 

Most methods for analysis of dynamical systems are somehow related to 
Lyapunov functions, and a Lyapunov function appears more or less explic- 
itly in most analysis conditions. Dissipativity analysis [212,81], absolute sta- 
bility [215, 164], and analysis based on integral quadratic constraints [139] 
can all be viewed as methods for Lyapunov function construction. Quadratic 
Lyapunov functions is the cornerstone in the analysis and design of robust 
linear controllers using semidefinite programming [27]. 

While efficient software for semidefinite programming has not appeared 
until quite recently, the simplex method for solving linear programming 
problems is more than 50 years old [47]. Researchers are since long aware 
of the benefits of deriving results that can be verified via linear program- 
ming. Computer algorithms for construction of piecewise linear Lyapunov 
functions appeared already in the late 70’s, and has continued to attract a lot 
of attention, see e.g., [33, 34, 140, 23, 142, 147, 61, 107, 93]. The focus has been 
on polytopic Lyapunov functions and uncertain linear systems. An impor- 
tant exception is the work [147] that considers poly topic Lyapunov functions 
for piecewise linear systems. Highly related to our approach is the stability 
analysis proposed in [118] in which piecewise linear Lyapunov functions 
(that may have affine terms and are not necessarily polytopic) were con- 
structed using so-called facet functions. 

When quadratic Lyapunov functions do not suffice, it is very natural to 
consider functions that are piecewise quadratic. For example, for the Lya- 
punov functions used in the Popov criterion are piecewise quadratic when 
the nonlinearity in the feedback connection is piecewise linear. The idea of 
“patching together” piecewise quadratic Lyapunov functions in the state 
space has also been used in the analysis of specific nonlinear systems, see 
[165]. To the best of our knowledge, this book (and the parallel work [158, 
71]) is the first that presents a systematic methodology for computation of 
piecewise quadratic Lyapunov functions for general piecewise linear sys- 
tems. 

Lyapunov Function Computations in the Equality-Constrained Format 

Rather than using the natural equality-constrained description of piecewise 
quadratic Lyapunov functions (4.13) we have introduced a special matrix 
format to enforce continuity of the Lyapunov function candidate. This for- 
mat enables the Lyapunov function search to be carried out using widely 
available software for semidefinite programming and does not force the user 
to rely on a solver that supports equality constraints. Although this issue 



82 




4.12 Comments and References 



was critical for the early developments of the results in this book (when no 
equality constrained solvers were publically available) it is related to op- 
timization technology that might change and, to some extent, already has 
changed. 

In the equality-constrained format, continuity of the Lyapunov function 
is enforced by introducing the equality constraints 

Pi = Pj + ejjtfj + tijefj Vj such that X* n Xj ± 0 (4.41) 

in variables Pj, Pj and Uj in the optimization problem. The ability to handle 
equality constraints also allows regions that contain the origin to be treated 
in an alternative way compared to Theorem 4.1. Let x eq be an equilibrium 
point in the interior of cell Xj. Then, the dynamics can be written as 

x = AiX a,i = Ai(x — x eq ) for XGlj 

and the Lyapunov function must be on the form 

(x - X eq ) T Pi(x - Xeq) + f j for X £ Xj 

for some f * e R. To accommodate for this in the Lyapunov function search, 
replace the continuity constraint (4.41) by 

Pi Pi%e q 

-xJ q Pj xJ q PjX e q + fj 

and replace the Lyapunov conditions 

Pi — E/T WiEi > 0 Aj Pi + PiAi + EfUiEi < 0 

by 

Pi — Ej WiEi > 0 Aj 1 Pi + PiAi + E. J WiEi < 0 

Here Ei is obtained from the matrix the matrix 

Gi | GiX eq + gi 

by eliminating every row whose last entry is non-zero (this is Algorithm A.l 
applied to the original cell identifiers of cell Xj). Note that the equality con- 
straints must be eliminated in order for the analysis inequalities to admit a 
solution with strict inequalities. By following the above procedure, all analy- 
sis results in this book can be converted to the equality-constrained format. 
More details on the equality-constrained format, as well as computational 
experience, can be found in [72, 71] (see also [172]). 



= Pj + efjtfj + Ujejj Vj such that Xj n Xj A 0 
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Sliding Mode Analysis and C 1 Lyapunov Functions 

A drawback of using Lyapunov functions that are continuous, but not neces- 
sarily continuously differentiable, manifests itself in systems with attractive 
modes. In order to ensure stability of Filippov solutions, we then have to an- 
alyze a possibly large differential inclusion for each sliding set. If we restrict 
the Lyapunov function to be C 1 this step is no longer necessary since 

dV(x,k) fA , dV(x,i) fA , WI 

— 7^ — {A&X + a^} — — — — {AkX + a&} \/k G K(i) 

The standard LMI conditions guarantee that the expression on the left-hand 
side is decreasing also on cell boundaries, so (4.28) is automatically satisfied. 

Since quadratic functions are C 1 the existence of a globally quadratic 
Lyapunov function guarantees exponential stability of Filippov solutions. In 
addition, it is straight-forward to restrict the piecewise quadratic Lyapunov 
function candidates to be C 1 both in the equality-constrained format and 
when using the special matrix parameterization introduced in Section 4.5. 
In the equality-constrained setting (4.13) is replaced by 

Pj = Pi I tij (hij + hjj) 

with Uj G R. For the matrix parameterization, without loss of generality, 
assume that the continuity matrices are on the form 



|_T n xn OnxlJ 

Then, by restricting the matrix parameter T to be on the form 



V 0 
0 T 0 



with T' diagonal and To symmetric, every function on the form (4.16) is C 1 
and piecewise quadratic. However, this restriction also limits the flexibil- 
ity of the Lyapunov function candidate. For example, this approach fails to 
establish exponential stability of the system studied in Example 4.10. 
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5 

Dissipativity Analysis 



A fundamental idea in systems and control is to view complex systems 
as the interconnection of simpler subsystems. Such a perspective is often 
helpful in bringing insight into and understanding about a dynamic sys- 
tem. When viewing a system as the interconnection of its components, it is 
natural to ask whether the analysis of a complex system can be based on 
the (hopefully simpler) analysis of its components. This is the idea behind 
input-output analysis, which has been a very successful tool in system the- 
ory. Roughly speaking, the idea is to replace detailed models of system com- 
ponents by input-output relationships (typically relating the energy of the 
input to the energy of the output) and then derive results for interconnec- 
tions of such models. The most well-known results may be the small-gain 
and the passivity theorems. Both allow stability of a feedback interconnec- 
tion to be verified from the analysis of its components. Hence, by establish- 
ing £ 2 gain or passivity properties of piecewise linear systems we can hope 
to establish stability of interconnections by invoking small gain and passiv- 
ity theorems. This chapter will provide such tools. 

An interesting aspect of this approach is that it allows us to use differ- 
ent tools for analyzing different components. For example, physical insight 
may help us to establish passivity of one subsystem and piecewise linear 
techniques can allow us to prove strict passivity of another subsystem. This 
allows us to analyze systems that combine piecewise linear subsystems with 
components that are not easily dealt with using piecewise linear techniques. 
Time delays is one example of such a component. 

As we have seen in Chapter 2, several important interconnections of 
piecewise linear systems are themselves piecewise linear. Thus, the tools 
developed in this chapter will give us the choice to either analyze a complex 
piecewise linear system as it stands, or to analyze the subsystems in isolation 
and then apply small-gain and passivity results. This allows us to trade-off 
complexity in the computations against conservatism in the analysis. 
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Chapter 5. Dissipativity Analysis 

5.1 Dissipativity Analysis via Convex Optimization 

Dissipativity is a very useful notion in the study of performance and robust- 
ness of dynamical systems. Roughly speaking, dissipativity means that the 
system absorbs more energy from its environment than what it supplies. In a 
more abstract setting a dissipative system is defined as a system that admits 
a supply rate (defining “input power”) and a storage function (measuring 
“stored energy”) so that energy is always dissipated [212, 81]. Many impor- 
tant system properties, such as £ 2 -induced norms and passivity, correspond 
to different supply rates. 

There is a close relationship between Lyapunov functions and storage 
functions. In some cases, the storage function will qualify as a Lyapunov 
function that can be used prove system stability. As exponential stability 
of a linear system implies the existence of a quadratic Lyapunov function, 
dissipativity of a linear system with respect to a quadratic supply rate im- 
plies the existence of a quadratic storage function [212]. With the develop- 
ments from the previous chapter at hand, it is natural to base a dissipativity 
analysis of piecewise linear systems on storage functions that are piecewise 
quadratic. Before giving some precise results we will illustrate the ideas on 
the problem of estimating the £ 2 -norm of a piecewise linear system. The 
initial step is based on a simple and transparent Lyapunov technique devel- 
oped in [212, 81, 199]. 

Performance Bounds from Dissipation Inequalities 

Consider the problem of estimating an upper bound on the £ 2 -induced gain 
from u to y of the system (2.3). In other words, we want to determine a 
constant 7 such that 

[ \y{t)\\dt < 7 2 f h(t)\\ Idt \/u(t) 

Jo Jo 

holds for all T > 0 . We will assume that x( 0 ) = 0 . The inequality can 
be verified if we can find a non-negative storage function V(x) >0 with 
y(x( 0 )) = 0 such that 

JjV(x{t)) < 7 2 u T u - y T y (5.1) 

along system trajectories. Integration of this inequality gives 

V(x(T)) -V(x(0)) < 7 2 / \\u(t)\%dt- f \\y(t)\\l dt. 

Jo Jo 



86 




5. 1 Dissipativity Analysis via Convex Optimization 



Since V(x(0)) = 0 and V(x) > 0, we have 

0<7 2 / h(t)\\ldt- f \y(t)\ 2 2 dt 
J o Jo 

and the desired bound follows. 

The central difficulty in applying the technique is to find a storage func- 
tion V (x) that satisfies the dissipation inequality (5.1). We will consider piece- 
wise linear systems and piecewise quadratic storage functions. This will al- 
low us to compute estimates on the £ 2 -gain using convex optimization. For 
a given partition a best upper bound can then be obtained by minimizing 7 2 
subject to the relevant inequalities. 

For linear systems, it is sufficient to consider quadratic storage functions 
and the exact £2 -gain can be found using the above procedure [27]. A re- 
sult of similar elegance for piecewise linear systems does not, to the best of 
the author’s knowledge, exist. In general, the techniques developed in this 
chapter will only return bounds on the system performance. 



Dual Bounds from Worst Case Disturbances 

Since we are working with bounds rather than exact solutions, it is useful to 
have measures on how good the computed bounds are. For the dissipativity- 
like computations in this chapter, such bounds can often be obtained by con- 
structing worst-case inputs. As the optimal value of the £ 2 -gain can be ob- 
tained from the solution of the Hamilton-Jacobi-Bellman equation 

dV(x) . . >. T T 

— {AiX + <H + BiU) = ju u-y y, 

ox 

it is natural to try to construct a worst case disturbance that attains equality 
in (5.1). Let V (x) be a solution to (5.1) and let 7 be the gain estimate obtained 
in this way. A lower bound on the £ 2 -induced gain can then be obtained by 
maximizing the expression 



dV(x) 

dx 



( AiX + at + Biu) + y T y - 7 2 u T u 



with respect to u . Simulating the piecewise linear system with this input 
and comparing the input and output norms then gives a lower bound on 
the £ 2 -induced gain. 



Successive Refinements via Upper and Lower Bounds 

The bounds obtained by piecewise quadratic computations often give sig- 
nificant improvements over computations based on globally quadratic stor- 
age functions. Moreover, by refining the state-space partition, it is possible to 
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introduce more flexibility in the piecewise quadratic storage functions. With 
increased flexibility the computations can be repeated in hope of achieving 
better estimates. As increased flexibility comes at the price of increased com- 
putations, the upper and lower bounds can serve as an aid in the trade-off 
between precision in the analysis and the cost of computations. 



5.2 Computation of £ 2 -induced Gain 



The analysis outlined above can be directly combined with the develop- 
ments in Chapter 4 to give LMI conditions for £2 -gain computations for 
piecewise linear systems. After verification of stability, for example using 
Theorem 4 . 1 , an upper bound for the gain can be obtained as follows. 



Theorem 5.1— Upper Bound on C 2 Gain 

Suppose there exist symmetric matrices T, Ui and W* such that Ui and W* 
have non-negative entries, while Pi = FjTFi and Pi = FjTFi satisfy 



PiAi + Af Pi + CjCi + EjUiEi PiBi 

BfPi -l 2 I 

'PiAi + Af Pi + CjCi + EjUiEi PiBi ' 

S?P* 



for i e Jo 



for i e h 



Then every trajectory x(t) with x( 0 ) = 0, / 0 °°(|Ml2 + \\^\\ 2 ) dt < 00 satisfies 

f'OO f'OG 

/ \\y\ldt<^ \\u\ldt. 

J 0 Jo 

The best upper bound on the £ 2 induced gain is achieved by minimizing 7 2 
subject to the constraints defined by the inequalities. □ 

Proof: See Section B. 3 . 



In analogy with the previous section, a lower bound on the gain can be 
computed by the construction of a worst case disturbance. For this purpose, 
we will consider disturbances on the form u = LiX obtained by maximizing 
the expression 

2x T Pi(AiX + Biu) + \Cix\l ~ l 2 \ u \l 

with respect to u . The precise formulas for the resulting feedback gains Li 
will be given in the next chapter and are omitted here. In the above expres- 
sion, Pi and 7 come from the upper bound computation. Simulating the sys- 
tem with this control law and comparing the input and output norms gives 
a lower bound on the C 2 gain. 
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Example 5.1— Analysis of a Saturated Control System 
Consider the control system shown in Figure 5.1. The output of the system 
Gi(s) is subject to a unit saturation. The closed loop dynamics is piecewise 
linear, with three cells induced by the saturation limits u = ± 1 . We set r = 0 




Figure 5.1 Saturated control system. 

and estimate the £ 2 -induced gain from the disturbance d to the output y. 
With the transfer functions 

g 2 § — |— y 

Gl( ' S ' 1 ~ 16s 2 + s + 2 G2 ( s ) “ 4s 2 + 3s + 12 

we obtain the results shown in Table 5.1. 

Table 5.1 Various upper bounds on C 2 gain. 



Method 


Gain Estimate 


Quadratic Lyapunov function 
Lure function with positive 77 
Lure without constraints on 77 
IQC for monotonic nonlinearities 
Piecewise quadratic Lyapunov function 


No solution found 

No solution found 

37.63 

5.62 

5.54 



Here “Lure function” means a Lyapunov function of the form V(x) = 
x T Px-\-rf f 0 Cx sat (s)ds and “IQC for monotonic nonlinearities” means again 
estimate computed based on [218] using the toolbox [138]. A lower bound 
on the £ 2 gain, computed for the linear region, is equal to 5.52. □ 

This example is somewhat contrived, but it illustrates the differences that 
can be obtained from the various approaches. Apart from being useful in 
analysis of disturbance rejection properties, computation of £ 2 -induced gain 
can be used for establishing robust stability in presence of norm-bound un- 
certainties. We will illustrate this in an example at the end of this chapter. 
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5.3 Estimation of Transient Energy 



In order to demonstrate the use of partition refinements we will apply the 
central idea to the estimation of the “transient energy” 




x t Qix dt 



for x(t) G X i7 i e I 



of a piecewise linear system. This can be seen as an alternative to simula- 
tion. The value of this integral depends on the initial value and we would 
like our estimate to also be a function of the initial value. In this way, one 
computation will allow us to estimate the integral for every initial state. We 
assume that Qi = Qj have the zero interpolation property, i.e. 



x 1 QiX = x 



-T 



Qi 0 
0 0 



x = x T QiX 



for i G io 



The desired estimates can be obtained from the following minor modifica 
tion of the stability analysis in Chapter 4. 

Theorem 5.2— Upper Bound on Transient Energy 
Let x(t) g U ieiXi with x(oo) = 0 be a trajectory of the system (2.3) with 
u = 0 for t > 0. Consider symmetric matrices T and U it such that Ui have 
non-negative entries, while Pi = FfTFi and Pi = FjTFi satisfy 

+ Fj? U iEi i G Io , 

+ UiEi i G 1 1 . 



< mf a;(0) T P io a;(0). □ 

Proof: See Section B.3. 



0 > PiAi + Aj Pi + Qi 
0 > PiAi + AJ Pi + Qi 

Then 




A lower bound can be obtained similary, by replacing Qi by — Qi in the 
analysis. A solution to the resulting inequalities then implies that 

rOO 

x(0) T Piox(0) < / x T Qixdt. 

Jo 

A best lower bound estimated in this way can be obtained by maximizing 
V(xo) subject to the relevant constraints. Although the computations are 
optimized for a specific initial value, the computed function V(x) bounds 
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the value of the integral for all initial values on the partition. Furthermore, 
if Qi > 0 for all i e I, a solution to the LMIs of Theorem 5.2 also satisfies the 
inequalities in Proposition 4.4. Thus, if the computed V(x) is positive on the 
partition it will also qualify as a Lyapunov function for the system. 

The following example applies the results on the problem of estimating 
the output energy of a piecewise linear system and illustrates the use of 
partition refinements to obtain better and better estimates. 



Example 5.2— Transient in Flower Example 

Consider again the piecewise linear system defined in Example 4.3 and in- 
troduce the output y = x i, i.e., let 



Ci* 



1 




Vie J. 



Figure 5.2 shows a simulated trajectory (left) and the corresponding output 
(right). Consider the problem of estimating the “output energy”, 

POO poo 

/ \y(t)\l dt= x(t) T Cj Cixit) dt 

Jo Jo 

from a given initial value. This can be done by direct application of The- 
orem 5.2 by letting Qi = CfCi . The output simulated from x(0) = [1 0] 

(shown in Figure 5.2) has total energy / 0 °° ||y(£)||f eft = 1.88, while the bounds 
obtained from Theorem 5.2 for the initial cell partition give 

pOO 

0.60 < / | j/ 1 j dt < 2.50. 

Jo 

A possible reason for the gap between the bounds is that the level curves 
of the cost function can not be well approximated by a piecewise quadratic 
function on the given partition. To improve the bounds, we introduce more 
flexibility in the storage function by repeatedly splitting every cell in two. 
This simple-minded refinement procedure, illustrated in Figure 5.3, is re- 
peated three times yielding the bounds shown in Table 5.2. 

Note that the bounds on the output energy optimized for the initial state 
(1,0) match closely over the the whole state space, giving good estimates of 
the output energy also for other initial states. □ 



5.4 Dissipative Systems with Quadratic Supply Rates 

The results of the previous sections can be generalized in a natural way to 
validation of more general dissipation inequalities. The same technique that 
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%2 




Figure 5.2 Trajectory of a simulation (left) and corresponding output (right). 

X 2 X2 




Figure 5.3 Upper (full) and lower (dashed) bounds on the storage function computed 
in Example 5.2. The bounds get increasingly tight when we move from 8 cells (left) to 
1 6 cells (right) . 



was used in £ 2 -induced gain computations can be applied to verification of 
dissipativity with respect to arbitrary quadratic supply functions. We give 
the following result. 

Theorem 5.3— Validation of Dissipation Inequalities 

Consider symmetric matrices T, Ui and Wi such that Ui and Wi have non- 

Table 5.2 Lower and upper bounds for output energy estimated via Theorem 5.2. 



Number of Cells 


Lower bound 


Upper bound 


4 


0.60 


2.50 


8 


1.33 


2.18 


16 


1.65 


1.98 


32 


1.78 


1.88 
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negative entries, while Pi = FjTFi and Pi = FjTFi satisfy 



0 < 


Ci D- 

0 I 


T 

M 


' Q D- 
0 I 


~ 


' PiAi + AfPi + EjUiEi PiB- 

BjPi 0 _ 


[0 < Pi - EjUiEi 









for i e Iq and 



[o< 

\ 


' Ci Di 

0 / 


T 

M 


~Ci 

0 


D- 

1 . 


~ 


PiAi + AJ Pi + EjUiEi PiBi 
BjPi 0 _ 


[0 < Pi - EjUiEi. 











for i G I\. Then every trajectory x(t) with x(0) = 0 and ||w(t )||2 dt < oo 
satisfies 



0 <f 


y 0) 


T 

M 


y( s ) 


Jo 


«(*)_ 




_u{s)_ 



Vt > 0 



□ 



Proof: Follows similarly to Theorem 5.1. 

Note that Theorem 5. 1 is a special case of this result where 



Theorem 5.3 can be used to establish induced gain and passivity properties 
of piecewise linear systems that can be used in analysis based on the small 
gain or passivity theorems. This type of results open up many possibilities 
as they allow freedom in whether to incorporate nonlinearities in the system 
description or to replace them by energy inequalities and use interconnec- 
tion results. The following example illustrates some of the ideas. 

Example 5.3— Robustness Analysis via the Small Gain Theorem 
Consider the system shown in Figure 5.4. This is a linear system with a dy- 
namic uncertainty A and a nonlinear spring <£>(•). The uncertainty block A is 
assumed to have induced £2 -gain less than one and the spring has the non- 
linear characteristic shown in Figure 5.4. Nonlinear spring arrangements of 
similar type can for example be found in engine control systems, see [120]. 
We set u = 0 and consider the system defined by a transfer matrix O(s) with 
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Figure 5.4 System with nonlinear spring characteristic and unmodeled dynamics 
(left) . Detailed spring characteristic (right) . 



state-space realization 
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dt 
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and outputs z d = x\,z s = x 2 . 

The tools derived so far give a large flexibility in how to analyze sta- 
bility of this system. A first crude approximation is to use a norm bounds 
on both the A-block and on the spring characteristic and then apply the 
small gain theorem. However, this approach fails for the suggested example 
since the norm of the linear system exceeds one (the Hoo norm is 5.73). 
Another approach is to consider the feedback interconnection of the linear 
system and the spring as a piecewise linear system and only treat the A 
block as uncertain. We may then start by using global sector bounds on the 
spring nonlinearity and estimate the £ 2 -gain from Wd to Zd (see Figure 2.1 1). 
This approach estimates the £ 2 -gain to 1.12 and stability can not be estab- 
lished. Using the piecewise linear sector bounds in Figure 2.11 (middle) , The- 
orem 5.1 verifies that the gain from Wd to Zd is less than 0.25 and closed loop 
stability follows. □ 

This example demonstrates how induced £ 2 -gain computations can be used 
for robustness analysis of piecewise linear uncertain systems, and how par- 
tition refinements are useful for obtaining sufficient accuracy in the analysis. 
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5.5 Comments and References 



Piecewise Linear Systems and Integral Quadratic Constraints 

The verification of dissipation inequalities could easily be extended to treat 
more general integral quadratic constraints (IQCs). For example, integral 
quadratic constraints with a frequency dependent weight rather than a con- 
stant matrix M could be verified using Theorem 5.3 by first introducing a 
state space realization of the weight and then include this dynamics in the 
system description. Once a certain integral quadratic constraint has been 
verified for a piecewise linear component it can be used in the general frame- 
work for IQC-based system analysis developed in [139]. A more straightfor- 
ward approach would be to perform an IQC-based analysis directly in the 
piecewise linear framework. The analysis of piecewise linear systems inter- 
connected with components described by IQCs could potentially be very 
useful. A good starting point for a stability analysis based on piecewise 
quadratic Lyapunov functions could be the Lyapunov approach outlined 
in [106, Section 1.7]. 

Nonlinear H 0 c Control 

Interpreted in time domain, the linear H 0 0 problem is concerned with atten- 
uation of the £ 2 -induced gain from disturbances to outputs. Nonlinear H 0 0 
refers to the extensions of these techniques to nonlinear systems, see [90]. 
In this context, £ 2 -induced gains are often estimated using finite difference 
discretization of the dissipation inequalities [91]. 

Incorporation of Hard Bounds 

In some cases we know that the disturbances, in addition to having bounded 
energy, satisfy hard amplitude constraints. It is straight-forward to modify 
Theorem 5.1 so as to account for this information. If we describe the input 
constraints via 

G^x ~L Hiu ^ 0 for i £ Iq G^x T~ HiU ^ 0 for i £ I 



the first inequality condition in Theorem 5.1 becomes 



'PiAi + Aj Pi + CiCf 


PXi 




~Ei O' 


T 

Ui 


'Et O' 


BfPi 


—7 T I_ 




Gi Hi 


Gi H 



and the second is analogous. 
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Controller Design 



So far, we have focused on deriving methods for stability and performance 
analysis of piecewise linear systems. At the very least such methods can be 
used in control design procedures that iterate between design and analysis 
steps. In this chapter, we will proceed by extending the piecewise quadratic 
techniques to design of feedback control laws. 

We will start out by showing how the problem of finding a globally lin- 
ear state feedback that stabilizes a given piecewise linear system can be 
cast as a convex optimization problem. This approach is elegant, but re- 
strictive in at least two ways: it is based on globally quadratic Lyapunov 
functions and it generates discontinuous feedback laws. When we try to al- 
leviate these shortcomings, however, we will find that it is not so easy to 
use piecewise quadratic Lyapunov functions as the basis for control design 
while keeping convexity in the computations. We will present an approach 
to optimal control of piecewise linear systems with piecewise quadratic loss 
functions. Feedback laws can then be derived from the solution to the asso- 
ciated Hamilton-Jacobi-Bellman equation. By considering Hamilton-Jacobi- 
Bellman inequalities rather than equations, we will show how convex opti- 
mization and piecewise quadratic functions can be used to compute bounds 
on the optimal performance. In particular, we will show how a lower bound 
on the cost that can be achieved using any control can be computed via 
semi-definite programming. The computations also suggest a natural piece- 
wise linear feedback law that tries to achieve this cost. The performance of 
this control law can be evaluated using the dual Hamilton-Jacobi-Bellman 
inequality, yielding an upper bound on the achievable cost. The dual in- 
equality can also be used for control design. However, the resulting ma- 
trix inequalities are bilinear and the associated design problem can only be 
solved locally. In an example, we show how the upper and lower bounds 
can be combined and a continuous piecewise linear state feedback with tight 
bounds on the guaranteed performance can be designed. 



M. Johansson: Piecewise Linear Control Systems, LNCIS 284, pp. 96-106, 2003. 
© Springer- Verlag Berlin Heidelberg 2003 




6 . 1 Quadratic Stabilization of Piecewise Linear Systems 

6.1 Quadratic Stabilization of Piecewise Linear Systems 

We have already discussed how piecewise linear systems can be conserva- 
tively analyzed as linear differential inclusions using quadratic Lyapunov 
functions. The same approach can also be used to formulate a number of 
central feedback control problems in terms of linear matrix inequalities [27, 
180]. However, this approach is very restrictive since it does not allow for 
bias terms in the system dynamics and it does not use any partition informa- 
tion in the analysis. In this section, we will describe an approach by Hassibi 
and Boyd [72] that alleviates both these shortcomings. 

Consider the piecewise linear system (2.3) under linear state feedback 

u = —Lx. 

This results in the closed loop system 

x(t) = ( Ai - BiL)x(t) + a; xeXi i el. (6.1) 

The quadratic stabilization problem ammounts to finding a feedback gain L 
so that the closed-loop system (6.1) admits a quadratic Lyapunov function 

V (x) = x T Px 

Assume that we use ellipsoidal cell boundings, i.e., that we for each cell Xi 
have determined matrices Si and Si such that 

|| Six + Si \\ 2 <1 Mx eXi 

Then the closed-loop system is quadratically stable if we can find a postive 
definite matrix P = P T > 0 and positive scalars Ui > 0 such that 



f 0 > (A - BiLfP + P(Ai - BX ) 






l 0> 


'(Ai - BXfP + P(Ai - BX) Pa- 


+ Ui 


1 

JEP 

i 

to 

Co 

c*. 

1 


aJP 0 _ 


-SiSi 1 - sf 



These conditions are not jointly convex in the variables L, P and u it but 
convexity can be recovered using the following result. 



Lemma 6.1 

There exists a solution P = P T > 0 and u > 0 to the matrix inequality 

\A T P + PA-uC T C PB — uC T D 
" * ( PB - uC t D) t u(I - D t D) 
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if and only if there exists a solution Q = Q T > 0 and v > 0 to 

\AQ + QA t - vBB t —vBD t + QC t ~ 

> {-vBD t YQC t ) t v(I — DD t ) 

Moreover, the soutions are related via Q = P _1 and v = u~ x . □ 

Proof: See [27, Section 5.1]). 

Using the above result, and introducing the variable Y via Y = LQ, we 
arrive at the following result [72]. 

Theorem 6.1— Quadratic Stabilization 

If there exists a positive definite matrix Q = Q T > 0, positive scalars vi > 0, 
and a matrix Y such that 

{ 0 > QAj + AiQ - Y T Bj -BiY i€ J 0 

QAf + A i Q-Y T Bl-B i Y-v i a i af QSj-v i a i sJ 1 \ 

{QSj -ViaisJ) T Vi(I - Sisf) 

Then, the feedback u = —Lx with L = YQ~ X renders the piecewise linear 
system (2.3) exponentially stable. □ 

The above result extends directly to the design of piecewise linear state feed- 
back laws on the form 

u = —LiX x g Xi. 

One then simply replaces the common matrix Y in the above conditions by 
separate variables Yi for each region. This feedback law is typically discon- 
tinuous and might introduce attractive sliding modes in the closed-loop sys- 
tem. However, according to the discussion in Section 4.12 the existence of a 
common quadratic Lyapunov function guarantees that all potential Fillipov 
solutions tend to zero exponentially. 

The reader who is familiar with the work on LMI-based analysis and de- 
sign of linear uncertain systems might recognize the techniques from analy- 
sis of norm-bound LDIs (see, e.g., [27, Section 5.1]). Indeed, all results devel- 
oped for norm-bound LDIs can be used for analysis and design of piecewise 
linear systems along the lines of the discussion above. 

6.2 Controller Synthesis based on Piecewise Quadratics 

The quadratic stabilization approach is elegant and gives a significant de- 
crease in conservatism compared to an approach that models the piecewise 
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linear elements as a polytopic uncertainty (compare Corollary 4.1). Still, 
with the powerful machinery developed in the previous chapters it is natu- 
ral to attempt to derive synthesis procedures based on piecewise quadratic 
Lyapunov functions. Unfortunately, there does not appear to be any simple 
analogue to the change-of- variables Q = P _1 when we base our analysis on 
piecewise quadratics and it does not seem to be possible to find a convex 
formulation of the design problem. 

In this section, we will investigate two approaches that allow us to de- 
sign controllers based on piecewise quadratic Lyapunov functions. We will 
focus on optimal control problems with piecewise quadratic costs. The first 
approach formulates the design problems as an optimization problem over 
bilinear matrix inequalities (BMIs). These problems are not convex and one 
cannot hope to find efficient algorithms that are guaranteed to find the glob- 
ally optimal solution (or even a feasible solution) for all problem instances. 
However, there are many heuristic algorithms for solving BMIs that appear 
to work reasonably well in practice. The second approach, first suggested 
in [169, 170], uses convex optimization to compute a lower bound on the 
objective that can be achieved by any control and suggests a natural piece- 
wise linear control law that tries to achieve this cost. The control law gives 
an upper bound on the cost. If the upper and lower bounds obtained in 
this way are not sufficiently tight, one can use this control law as a starting 
point for the BMI-based design procedure. When combined, the upper and 
lower bound computations allow us to design piecewise linear state feed- 
backs with guaranteed (and often quite tight) performance bounds. 

Optimal Control and Hamilton-Jacobi-Bellman Inequalities 

Consider the following general form of optimal control problem 



minimize f 0 °° L(x, u)dt 



subject to 



±{t) = f(x(t),u(t)) 
x(0) = Xq 



It is well known that solutions of this problem can be characterized in terms 
of the Hamilton — Jacobi — Bellman (H-J-B) equation 



0 = inf 

u 



dV 

—f(x,u) + L{x,u) 



(6.2) 



In fact, by integrating the inequality 



dV 

0 < —f(x,u) + L(x,u) 



Mx,u (6.3) 



99 




Chapter 6. Controller Design 
and assuming that x(oo) = 0 , we get 



V(x 0 ) - V(0) 



f°° dV 

—f(x,u)dt < 
n ox 



L(x : u)dt. 



Hence, every V that satisfies (6.3) gives a lower bound on the optimal value 
of the loss function. Similarly, for a given control law u = l(x), the inequality 



0 > ^-/(x, l(x)) + L(x, l(x)) 



Vx (6.4) 



guarantees that 



V(x 0 ) - V(0) > 




L(x , l(x)) dt 



and any V that satisfies (6.4) gives an upper bound on the cost. In particular, 
V has then decay rate given by —L(x, l(x)), which is typically negative, so V 
may serve as a Lyapunov function to prove that the control law is stabilizing. 

We will consider the case where / is piecewise linear and L is piecewise 
quadratic. The control objective is then to bring the system (2.3) from a given 
initial state x(0) = xo to rest at x(oo) = 0 while limiting the cost 

r-OO 

j(x 0 ,u)= (x{t) T Qg t) x{t) + u(t) T Rg t )u(t)) dt. (6.5) 

Jo 

Here i(t) is defined so that x(t) e X^ t y We will assume that 



Qi 



Qi 0 
0 0 



for i g Jo (6.6) 



and that Qi and Ri are positive definite. 

Upper Bounds via Bilinear Matrix Optimization 

Consider the piecewise linear system (2.3) under piecewise linear feedback 

u = — LiX — k = —LiX x G Xi, i G / (6.7) 

with li = 0 for i G Jo. Then, an upper bound on the associated loss function 
(6.5) can be computed via (6.4) in analogy with Theorem 5.2. 
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Theorem 6.2— Upper Bound on Optimal Cost 

Assume existence of symmetric matrices T and U if such that Ui have non- 
negative entries, while Pi = FjTFi and Pi = FjTFi satisfy 

q > (A* — BiLi) T Pi + Pi(Ai — BiLi ) + Qi + Ej UiEi 

[ Li 

q ^ (A-i — BiLi) T Pi + Pi(Ai — BiLi ) + Qi + Ej UiEi 
[ Li 

Then, every trajectory x(t) of (2.3) under the control 
U i£iX it x(oo) = 0 and x(0) = x 0 G X io satisfies 

J(xq,u) < mf XQp io x 0 (6.8) 

□ 



Li 

-R7 



i G Iq 



A' ,eJi 

(6.7) such that x(t) G 



Proof: Similar to Theorem 5.2 via the Schur complement (Proposition A. 2). 

When the Li s are fixed, the matrix inequalities are linear in Pi and an 
upper bound on the control performance can be obtained via semidefinite 
programming. The formulation can be used for controller design if we con- 
sider both Pi and Li as variables. However, in this case the inequalities con- 
tain bilinear forms PiBiLi. Optimization problems involving bilinear matrix 
inequalities are non-convex in general, but a number of heuristic methods 
for solving such problems have been proposed (see, e.g., [64, 66, 74, 196]). 
One of the simplest methods is the V-K iteration, which sequentially fixes 
one set of variables (say, the value function V) and optimizes over the oth- 
ers (typically, the controller gains) , then fixes the second set of variables and 
optimizes over the first. The procedure is repeated until no further perfor- 
mance improvements can be made. 

Unfortunately, the V-K iteration is not easily adopted to solving the BMIs 
in Theorem 6.2 since the objective (6.8) does not depend on Li. Instead, we 
will use the linearization method suggested in [74]. The basic idea of this 
approach is to linearize the design inequalities around an initial solution 
(Pi 7 Li), by letting Pi = P® + SP if Li = + 6 Li and ignoring higher-order 

terms. Thus, the expressions 

( [Ai — BiLi) T Pi + Pi(Ai — BiLi) 

in the inequalities of Theorem 6.2 are replaced by 

(Al) T P? + P z °Al + (AlfSPi + SPiAl - P® Bid Li - (SLifBfP? 
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where A ic = Ai — BiL®, Li is replaced by L^-\-SL it and the other entries of the 
design inequalities remain unchanged. This approximates the inequalities in 
Theorem 6.2 by linear matrix inequalities in S Pi and 5 Li. To guarantee that 
the linearization is valid, it is useful to add a trust region constraint such as 

m\\<a\\P?\\ 

Finally, rather than minimizing the right-hand side of (6.8) we minimize the 
quantity Xq (SPio)xo. 

An interesting feature of the BMI-based design approach is that it is easy 
to include additional constraints, such as continuity and boundedness of 
the control law (see, e.g., [171]). The following example illustrates the tech- 
niques. 

Example 6.1 

Consider the following simple model of an inverted pendulum 



x± = x 2 

x 2 = —0.1x2 + sin(xi) + u (6.9) 

We are interested in applying the proposed techniques to find a feedback 
control that brings the pendulum from rest at the stable equilibrium (7r,0) 
to the upright position (0, 0) while minimizing the criteria 

f'OO 

J(x, u) == / Ax\{t) + 4x2 (r) + w 2 (r) dr. 

Jt = o 

A piecewise linear model of the system (6.9) can be constructed by finding 
piecewise affine bounds on the system nonlinearity sin(xi). For the purpose 
of this example, we divide the interval [-4,4] into five segments and com- 
pute the bounds illustrated in Figure 6.1 (left). 

The equilibrium point (0, 0) is unstable, but the system is quadratically 
stabilizable and Theorem 6.1 returns the stabilizing control law 



u = —Lx = —1.20xi _ 1.15x2 



Theorem 6.2 estimates the performance of this control law to be 

J(xo, —Lx) < 90.68 

This control law can be used as starting point for the linearization method 
outlined above. To illustrate how additional constraints can be incorporated 
in the design procedure, we will require that the control law is continuous. 
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sm(xi) 





Figure 6.1 The left figure shows bounds on the system nonlinearity. The right figure 
shows the synthesized continuous and piecewise linear state feedback law. 



We use the format for continuous piecewise linear functions proposed in 
Section 4.10, i.e., we will let Li be on the form 

Li = L T Fi iel (6.10) 

Here, L is the parameter vector and Fi are the same continuity matrices that 
are used to parameterize the Lyapunov function candidate. A couple of Al- 
terations of the linearization procedure gives a continuous piecewise linear 
feedback shown in Figure 6.1 (right) with estimated performance 58.75. □ 



Lower Bounds via Convex Optimization 

The non-convexity of the control design problem in the general case is some- 
what discomforting. When the iterations cease to make progress we have no 
way of telling if the associated solution is close to optimal or not. In this sec- 
tion, we will use the alternative Hamilton-Jacobi-Bellman inequality (6.3) to 
derive a lower bound on the objective that can be achieved by any control. 
This bound can be computed via convex optimization. The computations 
also give a natural piecewise linear state feedback that tries to achieve this 
cost. The method can be used by itself or together with the design methods 
described above. The following result shows how the lower bound compu- 
tation can be cast as a semidefinite programming problem. 

Theorem 6.3— Lower Bound on Optimal Cost 

Assume existence of symmetric matrices T and U it such that Ui have non- 
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negative entries, while Pi = F?TFi and P* = FfTFi satisfy 



PiAi + Af Pi + Qi - Ef UiEi PiBi 
Bf P t Ri 

PiAi + AJPi + Qi — Ej UiEi PiBi 
BfPi Ri 



Then, every trajectory x(t ) of (2.3) with x(t) G X, x(oo) 
X io satisfies 



i e Io 



i g h 



0 and x(0) = xq g 



J(xo,w) > supxQP io xo 

T,Ui 



□ 



Proof: See Appendix A. 

Theorem 6.3 gives a lower bound on the minimal value of the cost func- 
tion J that can be achieved by any control. Upper bounds are obtained by 
studying specific control laws, such as the one obtained by the minimization 

min (^^f{x,u) + L(x,u)^j . (6.11) 

If the H-J-B equation (6.2) holds, then every minimizing control law is op- 
timal. However, if only the inequality (6.3) holds (for example as a result 
of solving the matrix inequalities in Theorem 6.3) then there is no guarantee 
that the control law minimizing (6.11) is even stabilizing. Still, the minimiza- 
tion problem can be used as the starting point for definition of control laws 
that will be used in our further analysis. 

Exact minimization of the expression (6.11) can be done analytically in 
analogy with ordinary linear quadratic control, using the notation 

Li = RT l BjP i} Li = R- 1 BjP i , (6.12) 

Ai — Xi BiL i: — Ai BiLi, 

Qi = Qi X PiB i R i 1 Bj Pi , Qi = Qi + PiB i R i ^Bf Pi. 

The minimizing control law can then be written as 

u(t) = — LiX x G Xi. 

If the control law is stabilizing, an upper bound on the optimal cost can be 
obtained from Theorem 6.2. 
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This control law is simple but may be discontinuous and give rise to at- 
tractive sliding modes (which are disregarded in Theorem 6.3). A procedure 
for constructing continuous feedback laws was suggested in [170]. In ad- 
dition, one can always use the value function computed in Theorem 6.3 to 
initialize the BMI approach. It is then easy to enforce continuity of the con- 
trol law using, for example, the techniques from Example 6.1. The following 
example illustrates the use of the lower bound computations. 

Example 6.2— Control of Inverted Pendulum (continued) 
Consider again the control of the inverted pendulum described in Exam- 
ple 6.1. A direct application of Theorem 6.3 gives the lower bound 

54.23 < J(x o, u) 

Thus, the contol law computed in Example 6.1 is close to optimal. It is easy 
to verify that the control law (6.12) obtained from the lower bound com- 
putation results in a stable closed-loop system with no attractive sliding 
modes. An evaluation of the performance of this control law using Theo- 
rem 6.2 gives the upper bound J(x 0 , u) < 59.64. This is slightly worse than 
what was achieved by the (computationally more intensive) BMI-approach. 
We conclude that the control law computed in Example 6.1 satisfies 

54.23 < J(x 0 ,u) < 58.75. 

Note that although the bounds are optimized for one particular initial value, 
the computed functions bound the performance for all initial values that 
result in trajectories that remain in the partition. The level surfaces of the 
upper and lower bounds on the value function are shown in Figure 6.2 (left). 
In addition, the computed control law is evaluated on the pendulum model 
(6.9) by simulation. The value of the loss function computed in this way is 
J(xq, u) = 54.49. □ 



6.3 Comments and References 



BMI-Based Controller Design 

The use of the matrix inequalities in Theorem 6.2 for control design is, to the 
best of the author’s knowledge, new. However, BMI-based control design 
procedures for piecewise linear systems have been suggested before. Slup- 
phaug [190] derived BMI conditions for quadratic stabilizability of discrete- 
time piecewise linear systems, Rodrigues et.al. have developed a number 
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Figure 6.2 The left figure shows lower (dashed) and upper (solid) bounds on the 
optimal costs. Although the bounds are optimized for the initial value xo = (7r, 0), 
the bounds match well over the full partition. The right figure shows a simulation of 
the computed control law. 



of BMI-based techniques for stabilization of piecewise linear systems with 
multiple equilibria, and Mulder and Kothare [143] used piecewise quadratic 
Lyapunov functions and BMI-optimization to design anti-windup control 
systems. The description of the BMI-based design procedure in this book is 
very short and the interested reader is urged to seek out the above references 
for more details and alternative formulations. 

Optimal Control and the Hamilton-Jacobi-Bellman Inequality 

Trajectory optimization and design of optimal feedback laws are impor- 
tant problems in science and engineering, and a wide variety of numerical 
methods for solving optimal control problems have been developed (see, 
e.g., [12, 37, 185]). We have considered an alternative approach based on 
Hamilton-Jacobi-Bellman inequalities and convex optimization. These ideas 
have been further developed in [168]. We also note that the H-J-B inequality 
has a long history in optimal control of systems with discrete states, see, e.g., 
[ 20 ]. 

Piecewise Linear Quadratic Control 

Work related to the term piecewise linear quadratic control has appeared 
in [11, 213]. In [11], linear quadratic control problems for piecewise linear 
systems were addressed by solving Riccati differential equations and the 
optimum had to be recomputed for each new final state. The reference [213] 
treats control design for linear systems with bounded controls. A piecewise 
linear control law is obtained by scheduling a set of state feedback con- 
trollers designed via linear quadratic theory. 



106 



7 

Selected Topics 



The analysis procedures developed so far can be extended in many useful 
ways. In this chapter, we present four specific extensions. 

Our first extension is to demonstrate how piecewise quadratic Lyapunov 
functions can be used to estimate regions of attractions. For systems with 
finite stability regions, the analysis conditions derived in Chapter 4 cannot 
be verified globally and the techniques cannot be immediately applied. We 
propose a simple solution to this problem that addresses the central issues 
of how to restrict the analysis region and optimize the size of the estimated 
stability domain. We observe that this approach yields significantly better 
estimates than recent methods based on the circle and Popov criteria. 

We then investigate how piecewise linear system approximations and 
piecewise quadratic Lyapunov functions can be used for rigorous analysis 
of smooth nonlinear systems. We show how approximation errors can be 
accounted for so that a successful stability analysis of the piecewise linear 
approximation also guarantees stability of the underlying smooth system. In 
addition, we establish a converse theorem stating that, in principle, when- 
ever the underlying system is exponentially stable we can compute a piece- 
wise quadratic Lyapunov function that proves it. 

It is often necessary to refine an initial partition in order to achieve suf- 
ficient accuracy in the piecewise linear approximation and adequate flex- 
ibility in the piecewise quadratic Lyapunov function candidate. It is then 
natural to ask how such partition refinements can be made efficiently and 
automatically. As a third extension we devise a simple method for automatic 
partition refinements. The algorithm uses linear programming duality in an 
attempt to improve flexibility and accuracy where it is needed the most. 

Our final extension is to show how fuzzy systems can be viewed as piece- 
wise linear differential inclusions and analyzed using piecewise quadratic 
Lyapunov functions. This allows an important class of fuzzy control sys- 
tems to be analyzed using substantially more powerful methods than was 
previously available. 
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Chapter 7. Selected Topics 

7.1 Estimation of Regions of Attraction 



In many cases, a local stability analysis of an equilibrium point cannot be ex- 
tended globally This occurs when the equilibrium point has a finite region 
of attraction (so that trajectories that start outside this region never converge 
to the equilibrium point) but can also happen for other reasons. The follow- 
ing example illustrates an interesting situation. 

Example 7.1— Lyapunov Analysis of Saturated System 
Consider again the double integrator under bounded linear feedback 



X 1 = X2 

X 2 = — sat(2xi + 3^2 ) 

We know from Example 4.1 that the origin is locally stable and that a (rel- 
atively small) region of attraction can be estimated using Proposition 4.1. 
Although an attempt to establish global stability using Theorem 4.1 fails, it 
is easy to verify that the piecewise quadratic Lyapunov function 

i i /*2 xi+3cc 2 

V(x u x 2 ) = -x\ + - J sat (z)dz (7.1) 

proves global asymptotic stability of the origin. The failure of Theorem 4.1 
can be understood by noting that the convergence is not exponential in the 
saturated regions and that the Lyapunov function (7.1) cannot be bounded 
in the sense of Lemma 4.1 (the Lyapunov function grows quadratically in 
the x\ -direction, but only linearly in the x 2 -direction). □ 

Even if the stability conditions developed in Chapter 4 fail to hold on the 
natural partition of the piecewise linear system one can still hope to find a 
larger stability domain than the one resulting from a purely linear analysis. 
A simple way to extend the local analysis is to try to verify the stability 
conditions on some ellipsoid £^(r) = {x \ x T Px < r 2 } C X. By gradually 
increasing the radius r of the ellipsoid, one increases the analysis domain, 
and hopefully also the resulting domain of attraction. Let 



Sa(t ) 



-P 0 

0 r 2 



and express £a(?) as 

£a ( r) = {x | x T S A (r)x > 0} 



(7.2) 



Stability analysis on X n £a(?) can then be carried out by adding the relax- 
ation term Ui§A(r) to the LMI conditions of Theorem 4.1. This leads to the 
following straight-forward but useful extension. 
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Proposition 7.1 

Consider positive scalars Ui and wi and symmetric matrices T, Ui and Wi 
such that Ui and Wi have nonnegative entries, while Pi = FjTFi and Pi = 
FfTFi satisfy 

f 0 > Aj Pi + PiAi + Ej 1 UiEi 

< T ® G Jo 

lO < Pi - Ef WiEi 

( 0 > AJ’ Pi + PiAi + Ej 1 UiEi + uiSa(x ) 

< iG/i 

lO < - EjWiEi - WiS A (r) 

Then, every level set 

Pa = {x | x T PiX < a x £ Xi , i G 1} 

such that C £ A (r) is a domain of attraction of the origin for the system 
(2.3) in the sense that all trajectories of (2.3) with u = 0 and x(0) e V a tend 
to the origin exponentially. □ 

Since there might be several solutions to the conditions of Proposition 7. 1 , it 
is natural to try to maximize the “size” of the computed level sets. However, 
determining the volume of a convex body (even if it is restricted to be a 
convex polyhedron) is computationally hard in general [69] and a formula 
that allows direct optimization of the volume of V a appears out of reach. For 
quadratic Lyapunov functions, V (x) = x T Px, one alternative to maximizing 
the volume of V a is to minimize the trace of P (see Section A. 2). Inspired by 
this criterion, we suggest to minimize the sum of traces of the matrices Pi 
that describe the Lyapunov function in each region. Hence, we propose to 
solve the convex optimization problem 

minimize Em ( 7 - 3 > 

iei 

subject to the conditions in Proposition 7.1. This is only a heuristic criterion, 
but it appears to work very well in practice. 

It makes intuitive sense to apply Proposition 7.1 iteratively. By restrict- 
ing the analysis region to some ellipsoid Sa(v) we can compute a piecewise 
quadratic Lyapunov function V(x) = x T P®x and an associated estimate of 
the region of attraction. This estimate often captures the shape of the actual 
region of attraction better than the initial ellipsoid and an improved esti- 
mate can be obtained by restricting the analysis region to a scaled version of 
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this domain. This can be done by replacing S A {r) with 



[0 r 2 J 4 

in Proposition 7.1 and sweeping r to find the largest value for which the 
revised analysis conditions admit a solution. 

To extract the best region of attraction that can be estimated using Propo- 
sition 7.1 we need to find the largest level set of V ( x ) contained in the anal- 
ysis region. A good estimate can often be obtained as follows: 

Proposition 7.2 

For every i e I, let a * be the largest a* such that 



&a(t) > Wi 



0 

0 



0 

Oii 




+ GfWiGi 



has a feasible solution W* y 0, wi > 0. Then 



T> a C S A {r) 



for all a < min ie j a*. □ 

Proof: For x e V at the inequality implies x T S A (r)x > 0, and hence x e S A (r). 

The following example demonstrates how Propositions 7.1 and 7.2 can 
be used to compute significantly better estimates of the region of attraction 
for the saturated system than what can be obtained from Proposition 4.1. 

Example 7.2— Extended RoA Estimation for Saturated System 
By restricting the analysis to a ball (. P = I) and gradually increasing the 
analysis radius r, exponential stability can be established for a very large ra- 
dius (in the order or r = 1E4). The level surfaces of the computed Lyapunov 
function are shown in Figure 7.1. Note that this is a drastic improvement 
over the linear analysis suggested in Example 4.1. 

We can understand the improved results by noting that even if the Lya- 
punov function (7.1) cannot be bounded by quadratic functions globally it 
can be bounded by quadratics on any compact domain. Numerical precision 
appears to limit how far the analysis domain can be extended. □ 

The choice of initial analysis ellipsoid is crucial for obtaining good results 
using Proposition 7.1. A simple choice is to use the ellipsoid obtained by 
the linear analysis of Proposition 4.1 but better analysis ellipsoids can be 
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Figure 7.1 Although the system trajectories do not have uniform global exponential 
decay in the saturated regions (left) , it is still possible to find a piecewise quadratic 
Lyapunov function that verifies stability for a very large set of initial values. 



obtained for particular classes of systems. In the rest of this section we will 
focus on linear systems with saturation. We will show how the resulting ap- 
proach gives significant improvements over alternative methods and yields 
estimates close to the true regions of attraction. 

Estimating Regions of Attraction for Saturated Linear Systems 

Consider a linear system under saturated linear feedback 

x{t) = Ax(t) + bsat(k T x) x(0) = xo (7.4) 

where sat(-) denotes the unit saturation. For simplicity we will focus on 
the case of a single saturation nonlinearity (the multiple saturation case can 
be dealt with using similar ideas). A piecewise linear description for this 
system has been worked out in Examples 2.1, 4.5 and 4.4. 

In Example 4.1, we illustrated how a simple estimate of the region of 
attraction can be computed using a linear analysis restricted to the non- 
saturated region. However, as shown in [160, 82], it is possible to obtain 
much better estimates at a slight increase in computations. By modelling the 
saturation as a locally sector-bounded element 

q/r < sat(g) < q for \q\ < r 

(see Figure 7.2) we obtain a polytopic model of the system dynamics 

x(t) G co {(A + bk T )x(t ), ( A + r~ 1 bk T )x(t)J 

This model is valid for all trajectories x(t) of (7.4) that satisfy 

\k T x{t)\ < r \/t > 0 
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\ RaM 




Figure 7.2 Modeling the saturation nonlinearity as a sector-bounded element (left) 
and the associated analysis region Ra ( r ) in the state space (right) . 



An estimate of the domain of attraction can be found by searching for a 
simultaneous Lyapunov function for the extreme systems of the polytopic 
model (similarly to Proposition 4.2). Since the model is only valid within 
the analysis region 

Ra{ r) = {x | — r < k T x < r} 

any level set of V(x) contained within Ra(c) is a domain of attraction of 
the system (7.4). The following result shows how to find the simultaneous 
Lyapunov function with largest level set (in terms of trace) fully contained 
in the analysis region (cf. [82, 160]) 

Proposition 7.3— Ellipsoidal RoA for Saturated Linear System 
A solution P = P T to the convex optimization problem 

minimize TrP 

subject to (A A bk T ) T P A P(A A bk T ) < 0 

(A A r~ 1 bk T ) T P A P(A A r^bk?) < 0 (7.5) 

r 2 k T ] 

> 0 

k P 

guarantees that V r = {x \ x T Px < 1} is a region of attraction for (7.4). □ 

To find the largest domain of attraction that can be estimated using the for- 
mulation above, we need to do a search over the parameter r. The result- 
ing estimate is often substantially better than what can be obtained using 
a purely linear analysis. Moreover, since the estimated region of attraction 
typically captures the rough shape of the true stability domain it is natural to 
use V r as analysis ellipsoid in the piecewise quadratic approach. Thus, we 
propose to first apply Proposition 7.3, use the optimal solution P to define 
the matrices Sa in (7.2), and then apply Proposition 7.1. 
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Finally, when we restrict the analysis domain to the ellipsoid Sa(t) we 
can also refine the cell boundings to reflect the relevant parts of the saturated 
regions. Since x e £a{t) implies that \k T x\ < r, we can replace the cell 
boundings suggested in Example 4.5 by 



Ei 



'-k T -1 




’ k T -l' 


k T r 


P,3 = 


— k T r 



Our proposed methodology for finding a piecewise quadratic estimate 
of the domain of attraction for the linear system with saturation (7.4) can 
now be summarized as follows: 

Algorithm 7.1— RoA Estimate for Saturated Linear System 

1 . Compute an initial estimate of the domain of attraction using the ro- 
bust approach described in Proposition 7.3. 

2. Use this quadratic estimate to define an initial analysis ellipsoid via 
(7.2) and to refine the cell boundings for the saturated regions. 

3. Compute a piecewise quadratic estimate by minimizing the objective 
function (7.3) subject to the conditions in Proposition 7.1. Extract the 
estimated region of attraction using Proposition 7.2. 

4. If desired, use the piecewise quadratic estimate to redefine the analysis 
region and go to step 3. 

□ 



The foilwing two examples illustrate the approach. 



Example 7.3— [160] 

Consider the saturated linear system (7.4) defined by 



A = 


" 0 


1 " 


b = 


0 


k = 


" 2 " 




1 


0 




-5 




1 



The domain of attraction estimated by the various methods are shown in 
Figure 7.3. Note that the piecewise quadratic Lyapunov function matches 
the actual region of attraction very closely and that dynamics in the satu- 
rated regions have unstable equilibrium points in =b (5 0 ). In fact, the level 
sets of the computed Lyapunov function are parabolic rather than ellipsoidal 
in the saturated regions (the matrices Pi and P 3 both have one negative 
eigenvalue). □ 
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Figure 7.3 Domain of attraction estimated by Proposition 7.3 (dashed) and Proposi- 
tion 7.1 (full). The shaded region is the true region of attraction (obtained via simula- 
tion) . 



Example 7.4 

As a second example, consider (7.4) with 



A = 


" 0 


1 


b = 


0 


k = 


’ 2 " 




0 


0.1 




-1 




3 



In this example, we apply the iterative procedure where the initial piecewise 
quadratic estimate is used to define the analysis region in the second pass. 
We find the piecewise quadratic estimate shown in Figure 7.4 (full line). 
Again, the result is significantly larger than what is produced by Proposi- 
tion 7.3 and matches the true region of attraction very well. □ 
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Figure 7.4 The domain of attraction can be improved substantially using the iterative 
procedure (full line) . 



7.2 Rigorous Analysis of Smooth Nonlinear Systems 



7.2 Rigorous Analysis of Smooth Nonlinear Systems 

In this section, we will study how the piecewise linear techniques devel- 
oped in Chapter 4 can be used to analyze smooth nonlinear systems. We 
will consider an approach where the smooth nonlinear system is approx- 
imated by a piecewise linear system and show how approximation errors 
can be accounted for in the analysis to yield rigorous results for the un- 
derlying smooth system. In addition, we will formulate a converse theorem 
which states that whenever the smooth system is exponentially stable, we 
can compute a piecewise quadratic Lyapunov function that proves it. 

Intuitively, if a piecewise linear approximation of a smooth nonlinear 
system is “accurate enough” and the approximation admits a piecewise 
quadratic Lyapunov function, this Lyapunov function should also be valid 
for the underlying smooth system. To transfer this intuition into a formal 
statement we must take the approximation errors between the piecewise lin- 
ear model and the smooth system in account. One way of doing so is to use 
piecewise linear differential inclusions (c.f. Section 4.7). Another alternative 
is to use a norm bound of the approximation error as follows. 



Theorem 7.1— Norm Bound Approximation Errors 

Let x(t) be a piecewise C 1 trajectory of the system x = f(x) and assume that 

\\f(x) - AiX - Oi 1 2 < e<||x |2 xeXi,ieI. 

If there exists numbers 7 * > 0 , symmetric matrices Ui and W* with non- 
negative entries, and a symmetric matrix T such that Pi = FjTFi and Pi = 
FfTFi satisfy 

—2tciil > Af Pi + PiAi + UiEi ( 7 . 6 ) 

EjWiEi < Pi < lt I (7.7) 

for i € Iq and 



— 2e^7^/ > AjPi + PiAi + Ej 1 UiEi (7.8) 

EjWiEi < Pi < 7i J (7.9) 

for i € h, then x(t) tends to zero exponentially. □ 

Proof: See Section B.5. 
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In the introductory examples, we used a particularly simple approach for 
constructing a piecewise linear approximation of a smooth nonlinear func- 
tion. A compact and convex subset X C R n is partitioned into simplicies by 
specifying a number of grid points z/& (vertices of the simplices). A piecewise 
linear approximation is then constructed by matching the smooth function 
at the grid points, f{vk) = Mvk + a i> and interpolating linearly between 
these points. The approximation error induced by this scheme depends on 
the smoothness of the function / and the size of the simplices (see [151] for 
a precise statement). Hence, if no solution to the matrix inequalities in The- 
orem 7.1 can be found, it makes sense to refine the partition by introducing 
new grid points (that decrease the size of the simplices) and to try again. 
Such subdivisions simultaneously reduce the approximation error and in- 
crease the flexibility of the Lyapunov function candidate. 

It is natural to ask how restrictive this approach is compared to a theorem 
based on arbitrary continuous Lyapunov functions. The answer is given by 
the following result which demonstrates that, in priciple, whenever a Lya- 
punov function exists there also exists a norm-bound piecewise linear ap- 
proximation and a piecewise quadratic Lyapunov function that satisfies the 
relevant matrix inequalities of Theorem 7.1. 

Theorem 7.2— A Converse Theorem 

Let / G C 1 (X) where X is a bounded invariant polytope for the system 
x = /(#). If the system is globally exponentially stable on X, then for every 
sufficiently refined simplex partition X with corresponding matrices Ei and 
Fi defined Proposition A.12 and with Lb satisfying A^i = there exists 
a solution Ui , U i: Wi , W* and T to the inequalities in Theorem 7.1. □ 

Proof: See Section B.5. 

The partition refinements decrease the approximation error and increase 
the flexibility of the Lyapunov function candidates, but introduce additional 
matrix inequalities in the analysis problem. For such an approach to be use- 
ful in practice, it is important to have support for partition refinements so 
that increased flexibility is introduced only where it is needed. The problem 
of automated partition refinements is studied next. 



116 




7.3 Automated Partition Refinements 



7.3 Automated Partition Refinements 

When analyzing a piecewise linear system with a given partition it is nat- 
ural to use the same partition for the Lyapunov function. However, this is 
not a definitive choice: some systems admit globally quadratic Lyapunov 
functions while other systems require repeated partition refinements before 
the Lyapunov function candidate becomes flexible enough. For example, 
the initial partition for the dynamics in Example 4.11 had two cells but the 
Lyapunov function needed a refined partition of four cells before a solu- 
tion to the analysis inequalities could be found. To verify stability of the 
smooth nonlinear system in Example 4.8 it was necessary to refine a (glob- 
ally) linear differential inclusion into a piecewise linear differential inclu- 
sion. In Chapter 5, we illustrated how tighter performance bounds could be 
obtained by refining the partition in order to increase the flexibility of the 
piecewise quadratic storage function. The approach for analysis of smooth 
nonlinear systems described in the previous section relies on repeated par- 
tition refinements to decrease approximation errors and increase flexibility 
of the Lyapunov function candidate. 

For simple systems, the partition refinements can often be made in an 
ad-hoc manner. For more complex systems, however, it is important to have 
tools that indicate where partition refinements are needed the most. In this 
section we will show how linear programming duality can be used for au- 
tomated partition refinements in the stability analysis based on piecewise 
linear Lyapunov functions. 

Introducing Flexibility Where Needed 

To illustrate the ideas, consider the problem of finding a piecewise linear 
Lyapunov function on a polytopic partition. Rather than solving the linear 
programming problem as it stands in Theorem 4.4, we consider the follow- 
ing slight modification 

minimize r 

subject to r > pf Aii/k i £ Jo, v k e Xi (7.10) 

r > pJAiiy k i e Iu Vk € Xi 

In this formulation, exponential stability according to Theorem 4.4 is ob- 
tained when r < 0. If the optimal value of the linear program is positive, 
no piecewise linear Lyapunov function exists on the current partition. It is 
then reasonable to find the cell which imposes the strongest constraint on 
the optimization problem and to subdivide this cell in order to increase the 
flexibility of the Lyapunov function candidate. The computations can then 
be repeated, proving stability or suggesting further partition refinements. 
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The sensitivity information about which constraints restrict the optimal 
value r can be obtained from the solution to the associated dual problem 
(see, e.g., [21]). Since there are many constraints (i.e., many dual variables) 
associated to each cell, we propose to compute the “total constraint cost” for 
each cell as the sum of the dual variables associated to it. The cell with the 
largest constraint cost can then be subdivided, and an iterative refinement 
procedure can proceed along the following steps. 

Algorithm 7.2— Automated Partition Refinements 

1. Solve the linear program from Theorem 4.4, modified as in (7.10). 

2. If r < 0, the procedure has terminated and asymptotic stability has 
been proven. Otherwise, refine the cell with the highest constraint cost 
and return to 1 . 

□ 



Cell Splitting 

After deciding which cell to subdivide, one must also decide how this subdi- 
vision should be done. This appears to be a delicate issue, since not every 
splitting operation increases the flexibility of the Lyapunov function. The 
simple idea of splitting a cell by introducing a new vertex in its center has 
the disadvantage that the boundaries of the cells are never refined (see Fig- 
ure 7.5a). We therefore suggest to split a cell by introducing a new vertex 




Figure 7.5 Procedures for subdividing a simplex; insertion of a new vertex in the 
center of a cell (left), or at the center of the longest edge. 



in the center of its longest edge (see Figure 7.5b). Note that this operation 
induces a subdivision also of the neighboring cells. Since the vector field in 
cells containing the origin is homogeneous, we propose to split these cells 
by introducing a new vertex in the largest edge of the facet that does not 
contain the origin. 

The following example illustrates the proposed approach. 
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Example 7.5 

Consider the following system 
(x\ X'2 

[;i'2 = —ox I ~ X'2 - sat(xi + X'2 ) 



(7*11) 



where sat(-) denotes the unit saturation. This system is piecewise linear and 
has oscillatory dynamics in both the linear and the saturated operating re- 
gions. We know from Chapter 4 that piecewise linear Lyapunov functions 
may require a rather fine partition of the state space in order to prove sta- 
bility of oscillatory systems* Indeed, for the coarse initial partition shown in 
Figure 7.6 (left), no piecewise linear Lyapunov function can be found* Based 
on this initial partition the automatic refinement procedure terminates with 
the partition shown in Figure 7.6 (right)* The corresponding Lyapunov func- 
tion, shown to the left in Figure 7.7, guarantees exponential decay of all tra- 
jectories within the estimated region of attraction, shown as the outermost 
level set in Figure 7.7 (right). 




Figure 7.6 Initial partition (left) and automatically refined partition (right). 




Figure 7.7 Lyapunov function (left) and guaranteed region attraction (right). 



□ 



119 






Chapter 7. Selected Topics 

7.4 Fuzzy Logic Systems 



While fuzzy control was quickly embraced by industry, its acceptance in 
academia has been slow. One reason for this is the evident lack of systematic 
methods for analysis and design of fuzzy control systems. In this section, we 
will demonstrate how an important class of fuzzy systems can be modeled 
and analyzed as piecewise linear differential inclusions, and show how the 
pwLDI description can be constructed systematically from a given fuzzy 
rule base. This allows an important class of fuzzy systems to be analyzed 
using substantially more powerful methods than was previously available. 

Takagi-Sugeno Fuzzy Systems 

A large amount of work in fuzzy systems literature has been devoted to 
analysis of Takagi-Sugeno systems [194]. The behavior of these systems are 
described by a set of rules on the form 

Ri : IF x\ is F i: i AND . . . AND x n is F* ?n 

THEN x = AiX + ai, i =3 1, . . . , L. 

Normally, the fuzzy sets (describing the kth input variable in the Ah 
rule) are labeled using linguistic terms such as “Small”, “Hot” or “Satu- 
rated”. Contrary to classical logic, where a proposition such as “ x k is F^&” 
can only take on values 0 or 1, fuzzy logic allows propositions to be fulfilled 
to any degree in the interval [0, 1]. To support this, one introduces member- 
ship functions 



IH,k(x k ) : M 1 — > [0, 1], 



that to each xu assigns a degree of validity of the proposition “x& is F^”. 
Using extensions of the inference and compositional rules of classical logic, 
fuzzy logic provides a framework for reasoning using linguistic rules and 
information described by fuzzy sets. In particular, given a set of linguistic 
rules on the form (7.12) and a numerical observation x, fuzzy logic can be 
used to infer a numerical value of x. This allows us to specify dynamic sys- 
tem models using a set of linguistic rules, the associated membership func- 
tions and some inference parameters (see [50, 209, 7] for further details). By 
the appropriate restrictions of the fuzzy inference parameters [194, 195, 209] 
the dynamics infered from the rules (7.12) can be written as 

L 

x = y^/^(x) {AiX-\- di} (7.13) 

i= 1 
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where gi{x) are normalized membership functions defined as 



= n Fi,k( x ) 



k= 1 



W (x) = 



(7.14) 



The normalization implies that fii(x) satisfy 0 < m(x) < 1, J2iFi( x ) = 1* 

Drawing upon the work on quadratic stability and quadratic stabiliza- 
tion, the formulation (7.13) has been used as a basis for solve analysis and 
control problems for fuzzy systems using LMI computations, see e.g., [219, 
195]. Although a major breakthrough in the analysis of fuzzy system, these 
methods view the fuzzy system (7.13) as a linear differential inclusion and 
are subject to the same criticism as the quadratic stability analysis of piece- 
wise linear systems. First, they disregard the partition information encoded 
in the membership functions and the premises of the fuzzy rules. Second, 
they do not admit affine terms in the consequent dynamics. This is a ma- 
jor shortcoming, since many applications use affine Takagi-Sugeno systems 
(see, e.g., [192, 10]) and the function approximation capabilities of Takagi- 
Sugeno systems are significantly improved when offset terms are allowed 
[53]. Finally, they use a very restricted Lyapunov function class (quadratics). 

It is thus natural to extend the piecewise quadratic analysis to fuzzy sys- 
tems. Such an extension would allow fuzzy systems with affine consequent 
dynamics to be analyzed using very powerful Lyapunov functions in a way 
that accounts for the structural information expressed by the linguistic rules. 



A Piecewise Linear Perspective 

In order to clarify the link between fuzzy systems and the piecewise linear 
systems considered in this book, it is fruitful to consider fuzzy systems as 
a particular instance of operating regime based models [94, 144]. Operat- 
ing regime based modeling is a common name for techniques that construct 
a globally valid model of the system dynamics by combining simple local 
models, each valid within a certain operating regime. In this context, the spe- 
cial feature of fuzzy systems is that prior knowledge of operating regimes 
and locally valid dynamics is encoded using fuzzy rules. Each rule premise 
defines an operating regime and the associated rule consequent specifies the 
local model valid within this region: 



Ri : 



IF #i is Fi t i AND . . . AND x n is 

operating regime specification 

THEN x = AiX + a>i, i = 1, . . . , L 

local dynamics 
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From the representation (7.13) we see that in regions where pi(x) = 1 for 
some i, all other normalized membership functions evaluate to zero and the 
dynamics of the system is given by x = Aix + a^. We will call this region 
the operating regime of model i. Between operating regimes there are regions 
where 0 < gi(x) < 1. In these regions, the system dynamics is given by 
a convex combination of several affine systems. We will call these regions 
interpolation regimes. As we will see next, the partitioning into operating and 
interpolation regimes is often polyhedral and can derived directly from the 
fuzzy rule base. 

The Geometry of Fuzzy Partitions 

There is more structure in fuzzy system partitions than what is directly vis- 
ible in the formulation (7.13). Consider for simplicity the case when the 
model scheduling is governed by one variable, say, x k : 

Ri : IF x k is Fi THEN x = AiX + a* i s= 1, . . . , L. (7.15) 

The operating regimes where /^(x*.) = 1 induce intervals in the schedul- 
ing space and parallelotopes in the state space. An example of membership 
functions and the associated partitioning is shown in Figure 7.8. 





Figure 7.8 The normalized membership functions (left) and the resulting partition of 
the state space into operating and interpolation regimes (right) . 

A similar partition structure results when rules are formed using the 
AND connective in a higher dimensional scheduling space. The rules then 
take the form (7.12), and the normalized membership functions are obtained 
as in (7.14). Since gi(x) = 1 in operating regime i, we must have pi t k(%k) = 1 
for all k = 1, ... , n. Hence, operating regime i is given by the intersection of 
the cells induced by the fuzzy sets that describe each propositional variable 

n 

Xi = {x\ m(x) = 1} = P| {x | lH,k( x k) = 1} 

m = 1 
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This allows us to view the partition resulting from rules formed with the 
AND connective as the composition of several simple partitions, see Fig- 
ure 7.9. The induced operating and interpolation regimes are also in this 
case parallelotopes and their description can be obtained directly from the 
membership function of the simple propositions. 



Figure 7.9 The fuzzy partition with scheduling in two variables (bottom) can be de- 
rived from the intersection of simple partitions induced by each propositional variable 
(top and center) . 



Fuzzy Systems as Piecewise Linear Differential Inclusions 

The discussion above establishes that affine Takagi-Sugeno systems can be 
modeled as piecewise linear differential inclusions. The fuzzy rules (7.12) in- 
duce a partitioning of the state space into convex polyhedral sets that act as 
operation or interpolation regimes. Within each cell, the dynamics is given 
by a state-dependent convex combination of affine dynamical systems 



o 

6 



0 

6 



0 

e 





k£K(i) 



where 0 < g k (x) < 1 , Y^keK(i) = h and 



K(i) = {k | /xfc(x) > 0 for some x G Xi\ . 
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In particular, if Xi is an operating regime, then K(i) contains only one el- 
ement. By disregarding the state dependence of the membership functions 
within each cell, we can embed these systems into the class of pwLDIs 

x G co {Akx\ x G Xi 

keK(i) 

and Theorem 4.2 applies directly. This allows an important class of fuzzy 
systems to be analyzed using much more powerful methods than was pre- 
viously available. 

If the associated LMI problem has too many variables, one might con- 
sider alternative formulations based on solving the analysis problem in two 
steps or the use of quadratic cell boundings (compare Section 4.9). Since the 
cells are often parallelotopes, quadratic cell boundings of minimal volume 
can be computed efficiently using Proposition A. 10. 

An Example 

In order to demonstrate the feasibility of the approach to problems of real- 
istic size, this section presents a piecewise quadratic stability analysis of a 
25-region fuzzy system. The system dynamics is given by the nine rules in 
Table 7.1. The membership functions of the fuzzy propositions “x* is Fi 



Table 7.1 Rule base for 25 region fuzzy system. 



IF 


x\ is negative 


AND 


^2 


is positive 


THEN 


X 


= A\x + a i 


IF 


x\ is zero 


AND 


^2 


is positive 


THEN 


X 


= A 2 x H- a 2 


IF 


x\ is positive 


AND 


a?2 


is positive 


THEN 


X 


= A 3 x + a 3 


IF 


x\ is negative 


AND 


a?2 


is zero 


THEN 


X 


II 

+ 

© 


IF 


x\ is zero 


AND 


a?2 


is zero 


THEN 


X 


= A$x + a 5 


IF 


x\ is positive 


AND 


a?2 


is zero 


THEN 


X 


= Aqx + a,Q 


IF 


x\ is negative 


AND 


a?2 


is negative 


THEN 


X 


= A^x + <27 


IF 


x\ is zero 


AND 


^2 


is negative 


THEN 


X 


= A s x + a 8 


IF 


x\ is positive 


AND 


^2 


is negative 


THEN 


X 


= Agx + ag 



are trapezoidal as shown in Figure 7.8, and the resulting partitioning of the 
state space into operating and interpolation regimes is shown in Figure 7.10. 
Note that the dynamics have bias terms in all regions that do not contain the 
origin, and that the system matrices associated to some of the operating re- 
gions are non-Hurwitz (these operating regimes, given by the two last rules 
of the rule base, are slightly shaded in Figure 7.10). Hence, the standard con- 
ditions for quadratic stability cannot be applied. 
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Figure 7.10 Partition of the fuzzy system defined by the rules in Table 7.1 into operat- 
ing regimes and interpolation regimes. The matrices Ai defining the dynamics in each 
operating regime are also shown. 




Figure 7.11 System response from a typical initial condition (left), and level surfaces 
of the computed Lyapunov function (right) . 
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As shown in Figure 7. 11 (left), simulations reveal a highly nonlinear be- 
havior but suggest that the system is stable. The linear matrix inequalities 
of Theorem 4.2 have a feasible solution proving exponential stability of the 
origin. The level surfaces of the computed Lyapunov function are indicated 
by dashed lines in Figure 7. 11 (right). We conclude asymptotic stability of 
the origin with a guaranteed region of attraction indicated by the outermost 
level set shown in Figure 7.1 1. 



7.5 Comments and References 

Estimating Regions of Attraction 

The problem of determining regions of attraction for nonlinear systems has 
a long history and a wealth of techniques have been developed in the liter- 
ature (see, for example, the survey papers [63, 24]). The problem is not only 
theoretically challenging, but also very relevant in practice. In many safety- 
critical systems (such as power systems, aerospace applications, and process 
control systems [152, 17, 179, 1]) it is essential to be able to determine the re- 
gion of safe operation. In addition, many high-performance controllers (e.g., 
[137, 217]) use time-optimal control for state transfer and a dynamic con- 
troller for regulation around a set-point. To be able to do safe switching be- 
tween the controller modes, it is important to know when the system state 
has entered the region of attraction of the dynamic controller (see [136] for 
more details). Although our approach applies to general piecewise linear 
systems, we have focused on linear systems with saturation for which both 
computational techniques and theoretical tools are more complete (see, e.g., 
[82, 160, 85, 174, 148, 65, 178, 86]). Techniques closely related to ours can be 
found in [82, 160, 85]. 

Establishing Convergence to a Set 

Some systems have equilibrium sets (multiple equilibrium points, limit cy- 
cles, etc.) rather than a single equilibrium point. Although the results devel- 
oped so far do not apply to this situation it is often possible to prove con- 
vergence to a set using a dual approach to the technique used for estimating 
regions of attraction. Let £a(t) be an ellipsoid. If we can find a Lyapunov 
function V(x) that satisfies the Lyapunov inequalities on X\£a(?) (i.e., on 
the partition with £a{t) removed) then any trajectory starting within the 
largest level set of V (x) contained in X converges to the smallest level set 
ofV(x) that contains £a{t). If £a{t) contains the origin, the Lyapunov func- 
tion search leads to similar matrix inequalities as in Proposition 7.1 but with 
Sa{t) replaced by —S^x)- 
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Stability Analysis of Fuzzy Systems 

The fuzzy community may have been first in using LMI computations for 
design and analysis of nonlinear control laws [219, 195] and the need to 
move beyond quadratic Lyapunov functions was recognized early [118]. In 
fact, the use of discontinuous piecewise quadratic Lyapunov functions has 
also been used in the independent work [54]. In comparison with these re- 
sults, an application of Theorem 4.2 to fuzzy systems constitutes a number 
of significant improvements. This includes the use of the 5-procedure to ex- 
ploit structural information and the ability to handle affine Takagi-Sugeno 
systems. Moreover, although [54] give procedures for computing different 
Lyapunov function expressions for each region, one must subsequently ver- 
ify certain boundary conditions to guarantee system stability. It is suggested 
that these boundary conditions could be verified by simulation. The ap- 
proach taken in this chapter avoids this non trivial step by parameterizing 
the Lyapunov function candidate to be continuous across cell boundaries. 

It should be noted, however, that when the normalized membership 
functions have global support (i.e., m{x) > 0 for all x e X) there is no clear 
partition of the state-space into operating and interpolation regimes, and 
the procedure described in Section 7.4 results in operating regimes where all 
consequent dynamics are active. An alternative approach for dealing with 
this class of systems has been proposed by Feng and Harris [55]. Inspired by 
the use of operating and interpolation regimes we can, for example, define 
cell Xi as 



Xi — {x | ^ ftlyk ip'k) yiy 

k 

Although rule Ri dominates in X ir all rules are active. Rather than using a 
differential inclusion that involves all consequent dynamics, Feng and Har- 
ris propose to re-write the dynamics as 



L L 

x = ^2 im(x) {AiX + <ii} = AiX + eg + ^ fJij(x) {A AijX + Aay} x e Xi 
i=l j=l 



where A A^ = Aj - Ai and A = a,j — a it Now, the derivative of the 
Lyapunov function V(x) = x T PiX along trajectories of the fuzzy system in 
cell Xi can be estimated as 

d ‘ L 



^ V (x(t)) = 2 x T Pi Ai + AK X ) AA; 



3 = 1 



< x T J Aj Pi + PiAi + 2 'y ' CjPi 

3 = 1 



(7.16) 
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where Qj = ||AA ^||2 The scalars Qj are readily computed, 

and the search for a piecewise quadratic Lyapunov function whose deriva- 
tive is negative can be formulated using the techniques from Chapter 4. 

Finally, we have only considered fuzzy systems of the Takagi-Sugeno 
type. Stability conditions for Mamdani-type fuzzy systems based on piece- 
wise quadratic Lyapunov functions can be found in [119]. 
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Linear Hybrid 
Dynamical Systems 



Hybrid dynamical systems are systems that combine continuous dynamics 
and discrete events. In practice, hybrid systems arise when some parts of a 
physical system can be modeled by differential equations while other parts 
are more conveniently described by logic or discrete-event systems. A sim- 
ple example is a continuous process controlled by switching logic (imple- 
mented by say a programmable logic controller). Although the piecewise 
linear systems considered so far can be regarded as hybrid systems (the cell 
index can be seen as a discrete state variable whose value changes when the 
continuous state crosses a cell boundary) the discrete state has a very pas- 
sive role in this case, and this model class exhibits few behaviors that cannot 
be understood using standard nonlinear systems theory. 

In this chapter we will consider linear hybrid dynamical systems, a class 
of piecewise linear systems with a more prominent hybrid nature. In con- 
trast to the piecewise linear systems studied so far, the discrete state is no 
longer directly determined by the continuous state. The discrete dynamics 
is described by a finite automoton whose state changes when the continu- 
ous state vector enters certain transition sets. Hence, the transition rules for 
the discrete dynamics depend on both the continuous and the discrete state. 
This model is simple enough to allow us to develop computational analy- 
sis tools yet rich enough to capture the essentials of many practical hybrid 
systems. In relation to the piecewise linear systems we can view this hybrid 
extension as piecewise linear systems with overlapping cells in R n . 

We will describe two computational approaches for stability analysis 
of linear hybrid dynamical systems. One is based on piecewise quadratic 
Lyapunov functions that are discontinuous but decrease in value at every 
switching instant, while the other approach uses Lyapunov functionals. 
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Chapter 8. Linear Hybrid Dynamical Systems 

8.1 Linear Hybrid Dynamical Systems 

A linear hybrid dynamical system can be represented as 

(x(t) = = A i{t) x(t) + a i{t) 

{ ( 8 . 1 ) 

[i(t) = )) 

Here x(t) e M n is the continuous state vector and i(t) el C Z is the discrete 
state variable. The differential equations describe the continuous dynam- 
ics while the algebraic equation models the discrete dynamics. The model 
associates one affine vector field for the continuous state to each value of 
the discrete state. The discrete dynamics is described by a finite automoton, 
whose state changes when the continuous state vector enters certain transi- 
tion sets. More precisely, a transition from the discrete state j to the discrete 
state k occurs when the continuous state x(t) enters the transition surface 

Sjk = {x | jj k x = 0} (8.2) 

In this way the discrete state variable i(t) becomes a piecewise constant 
function of time, and the notation t~ in (8.1) indicates that i(t) is piecewise 
continuous from the right. A trajectory of the system (8.1) on [to,tf] is a pair 
(x(t),i(t)), where i(t) is piecewise constant and x(t) is absolutely continu- 
ous, that satisfies the model equations for almost all t e [to,tf]. 

In our analysis procedures we need to keep track of the feasible transi- 
tions in the discrete state. We represent each feasible transition j — > k by an 
ordered pair (j, k) and denote the set of all feasible transitions T. In addi- 
tion we allow affine mode invariants Ej to be specified for each discrete state. 
These play the role of cell boundings in the piecewise linear systems model 
and it is assumed that Ejx(t) >z 0 whenever i(t) = j. Similarly to the piece- 
wise linear case, the mode invariants do not influence the dynamics but are 
introduced to reduce conservatism of the analysis. The mode invariants can 
be derived from the transition surfaces and other invariance considerations. 

Linear hybrid dynamical systems can be conveniently visualized as di- 
rected graphs which we will call transition diagrams, see Figure 8.1. Nodes in 
the transition diagram correspond to discrete states and arcs describe pos- 
sible transitions between discrete states. There is an affine dynamics and a 
mode invariant associated to each node and a transition surface associated 
to each arc. We identify nodes by the corresponding value of the discrete 
state and represent arcs as ordered pairs (j, k) of nodes. The leftmost arrow 
in Figure 8.1 indicates the initial state. 

The following example illustrates the model class and the need to extend 
the piecewise quadratic analysis developed in Chapter 4. 
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= {% I fri x “ °} 




Figure 8. 1 Transition diagram of linear hybrid dynamical system in Example 8. 1 . 





Figure 8.2 Sample trajectory of the hybrid system projected onto the continuous state 
space (left) and corresponding time responses (right) . 



Example 8.1 

Figure 8.2 (left) shows a simulation of the system x(t) = A^xit) 



i(t) = 




if i(t ) = 1 and /p 2 #(t) = 0 
if i(t~) ■»= 2 and /^(t) = 0 



with i(0) = 1, switching boundaries 

n T 



f 12 — 



-10 -1 



and system matrices 
Ai = 



-1 -100 

10 -1 



/21 — 



2 




A2 



1 10 

-100 1 



(8.3) 



The simulations shown in Figure 8.2 indicate that the system is asymptot- 
ically stable. From the simulated trajectory it is also clear that the system 
cannot be analyzed using a Fyapunov function that disregards the influence 
of the discrete state. □ 

The state space of the model (8. 1) , R n x Z, can be thought of as a set of enu- 
merated copies of R n . From this perspective, a transition in the discrete state 
can be seen as the transfer from one copy of R n to the other, see Figure 8.3. 
The simulation in Figure 8.2 is the projection of this trajectory onto R n . 
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Figure 8.3 State space of hybrid system illustrated as a number of enumerated copies 
of R n . Changes in the discrete state transfers the state from one copy to the other. 



8.2 Analysis Using Discontinuous Lyapunov Functions 

The analysis procedures that we have developed for piecewise linear sys- 
tems have all been based on continuous Lyapunov functions. When we try 
to analyze the more complex linear hybrid dynamical systems, however, ad- 
ditional flexibility is often needed. One way of increasing the analytic flex- 
ibility is to consider Lyapunov functions that are discontinuous, but whose 
value decreases every time there is a change in the discrete state. As the fol- 
lowing simple argument shows, such Lyapunov functions can also be found 
by semidefinite programming. 

Let the Lyapunov function candidate be V(x,i) = x T PiX and assume 
that the discrete state initially have the value j. Consider the transition from 
discrete state j to k, which occurs when the continuous state enters the tran- 
sition surface Sj k . The requirement that the Lyapunov function be decreas- 
ing at the discrete state change can be expressed as 

x T PjX > x T P k x for x G Sjk = {x | fj k x = 0} 

which is equivalent to the following LMI in Pj,P k and ij k 
Pj ~ Pk + fjktjk + tjkfjk ^ 0. 

To formalize these ideas, we consider the problem of establishing expo- 
nential stability of the equilibium point x = 0. Let io C / be the non-empty 
set of discrete states for which x(t) = 0 is admissible, and let I\ = I\Iq. For 
each i g /, assume that the mode invariants are on the form 

EiX ^0 i G Iq? .EjX 0 i G L 
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To obtain strict inequalities in our analysis computations, we define 



fjk — 



I 0 



fjk 



if 0 G Sjk 
otherwise 



(8.4) 



and let Pi = [I 0] T Pi [I 0] for i e /q. We then have the following result. 



Theorem 8.1 

Consider symmetric matrices Ui and Wi with non-negative entries, symmet- 
ric matrices Pi and P ir and vectors tjk and tj k such that 



JO > Aj Pi + P%Ai + Ej UiEi 
\o <Pi-EjWiEi 


i e Jo 


(8.5) 


f 0 > Af Pi + PiAi + Ef UiEi 


jo KPi-EjWiEi 


i € h 


(8.6) 


0 < Pj ~ Pk + fjktjk + ^jk fjk 


(j, k) € T, j € h or k e h 


(8.7) 


0 < Pj ~ Pk + fjktjk + tjk fjk 


(j, k) e T, j, k€lo 


(8.8) 



then x(t) tends to zero exponentially for every trajectory of (8.1). □ 

There is a strong relation between Theorem 8.1 and Theorem 4.1. By allow- 
ing non-strict inequalities in (8.7) and (8.8), Theorem 4.1 can be seen as a 
special case of Theorem 8.1 where fj k = fji Vi, j G /. However, a formula- 
tion with non-strict inequalities is numerically very sensitive and most LMI 
solvers can not treat non-strict inequalities as they stand. Inherent algebraic 
constraints must first be eliminated. Theorem 4.1 can be seen as the outcome 
of such an elimination. 

Theorem 8.1 can be applied directly to the system of Example 8.1. 



Example 8.2 

Consider again the switching system (8.3). To illustrate the use of mode in- 
variants, we let 



Ex 



0 

0 



0 




-1 

-1 



The LMI conditions of Theorem 8.1 have a feasible solution 



Pi = 



17.9 

-0.89 



-0.89 

179 



P2 = 



739 

-38.1 



-38.1 

91.8 
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X2 V{t) 




Figure 8.4 Sample trajectory of the hybrid system (left) and the corresponding value 
of the Lyapunov function (right) . 



A simulated trajectory of the system and the corresponding value of the 
computed Lyapunov function are shown in Figure 8.4. The discontinuities 
in the Lyapunov function concur with changes in the discrete state. □ 



Extensions 

Several useful extension can be made to the above results. It is, for example, 
straight-forward to allow for nonlinear or uncertain dynamics by conisder- 
ing differential inclusions in each mode. Similarly to before, this results in 
multiple decreasing conditions (one for each extreme system defining the 
dynamics in each mode). 

It is also easy to extend the approach from considering global transition 
surfaces to allowing for polyhedral transition sets 

Sjk = {x I fj k x = 0 A G jk x y 0} 

In this case, (8.7) is replaced by 

0 <Pj —Pk + fjktjk + tjkfjk ~ GjkZjk&jk 

where Zj k has non-negative entries. The other inequalities follow similarly. 
The same approach can also be used for modeling uncertain transition sur- 
faces, since a solution to the inequality 

0 < Pj — Pk — Gj k Zj k Gjk 

guarantees that V(x,j) > V(x, k) for all x with G jk x >z 0 (see, e.g., [158]). 

Finally, it is possible to consider systems that allow for jumps in the con- 
tinuous state at discrete transitions. In particular, assume that we have affine 
jump-maps 



x(t) = H i{t - )i{t) x(t ) 



W 
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with Hu = 0. Then, the condition x T (t)Pjx(t) < x T (t~)Pix(t~) can be ex- 
pressed as 

0 < Pi — Hj k PjHjk + fjktjk + tjkfjk 

and the other inequalities follow similarly (note, however, that we need to 
require that the jump maps are linear when i, j e Io in order to have strict 
inequalities in our formulation). 



8.3 Stability Analysis Using Lyapunov Functionals 

An alternative approach to analysis of linear hybrid dynamical systems has 
been suggested by Hassibi, Boyd and How [73]. The approach is much more 
involved than the method presented in Section 8.2 and we will not present it 
in its fullest generality or detail here. Rather, we aim at giving the flavour of 
the approach, illustrate how it differs from the analysis based on piecewise 
quadratic Lyapunov functions, and demonstrate how it can be applied to 
the analysis of a simple hybrid system. More details can be found in [73, 71]. 

Lagrange Stability and a Lyapunov Functional Class 

Contrary to the analysis procedures derived so far, which consider exponen- 
tial stability, the analysis procedure described in this section is concerned 
with Lagrange stability. Recall that a dynamical system with state x is called 
Lagrange stable if \\x(t)\\ remains bounded for any given initial state x(0) — x 0 
[123]. The stability analysis conditions are based on the following result. 

Proposition 8.1— [71] 

Suppose that the functional V : (x, i) i— > R is such that H(x, i) —> oo as 
||x|| — ► oo. Suppose furthermore that V is bounded below, i.e., that there 
exists f3 > 0 such that V > —(3. Along trajectories of the dynamical system 
V(t) = V(x(t), i(t)) becomes a function of time t. If V(t) is a decreasing 
function of t, then the system is Lagrange stable. □ 

To make constructive use of Proposition 8.1, Hassbi et. al., propose a flex- 
ible class of Lyapunov functionals that can be computed via convex opti- 
mization. To simplify the presentation, we will describe the approach for a 
systems whose transition surfaces are parallell to each other 

Sjk = {x I fx - Qjk = 0} (8.9) 

In this case, the stability analysis is based on the Lyapunov functional 

V(t) = x(t) T Pi(t)x(t) + 2J (a(s)f T x(s) + [3(s))f T ^x(s)ds (8.10) 
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Here, the matrices P* describe a continuous piecewise quadratic Lyapunov 
function (like the ones we have used in Chapter 4). The parameters of the 
integral term, a(t) and /3(t), are piecewise constant and change values when 
a discrete state transition takes place, i.e., 



and 



= 



(3(t) = 



a jk 


if *(*-)=. 


a(t~) 


otherwise 


Pjk 


if »(*“) = . 


P(t~) 


otherwise 



In other words, there is one set of parameters (ajk,fljk) associated to each 
arc of the transition diagram. Note that the value of V(t) is not uniquely 
given by x(t) and i(t) but depends on the full history of the trajectory This 
contrasts the approach presented in Section 8.2, where V(t) is an explicit 
function of (x(t),i(t)). 

Computing Lyapunov Functionals via Convex Optimization 

We will now show how the conditions in Proposition 8.1 can be converted 
into constraints on the parameters of the Lyapunov functional (8.10) . To sim- 
plify the presentation even further, we will assume that the Lyapunov func- 
tion part is globally quadratic, i.e., l\ = [I 0] T P [I 0] for all modes. 

It is convenient to the Lyapunov function part be positive definite. To 
guarantee that the Lyapunov functional (8.10) is bounded below, it is then 
sufficient to require that 

1. The integral term remains bounded from below for each fixed i(t) = i. 

2. The net ammount of the integral term around any cycle of the transi- 
tion diagram is non-negative. 

To deal with the first requirement, assume that there is a transition from 
mode j to mode k at time £&. Then, for every t > tu such that the system 
remains in mode k we have 

(a(s)f T x(s) + (3(s))f T ^-x(s) ds = 

a jk (x(t) T ff T x(t) - g 2 jk )/ 2 + (3 jk (f T x{t) - g jk ) 

If|/ T x\ can become unbounded in mode k, we must require that > 0. 
To this end, we let I u C / be the set of modes for which \f T x\ can become 



136 




8.3 Stability Analysis using Lyapunov Functionals 



unbouded and require that 



a jk > 0 "ike I u ( 8 . 11 ) 

To satisfy the second requirement, suppose that x(t) hits the transition 
surface Sj k at time t k and subsequently hits S k i at time ti. Then, 

/ (a(s)f T x(s) -\-(3(s))f T — x(s) ds = / (a jk a + (3 jk ) da 
Jt k ds J 9jk 

= a jk(9 k l — 9j k )/2 + Pjk(9kl ~ 9jk) 

which is a linear expression in the parameters aj k and (3j k . For any cycle 
in the transition diagram, we can evaluate the net ammount of the integral 
term similarly. Let (Aq, & 2 , . . . , Aq) be a sequence of discrete modes that 
constitute a cycle in the transition diagram. Then, the condition that the am- 
mount of integral around the cycle be non-negative reads 



E 

i= 1 



a kik. 



Q kiJ r \ki~s r 2 Hkiki-^-i 



■i+i 



+ Pkik i+ 1 (0*4+1 * 4+2 0*4*4+ 1 ) > 0 (8.12) 



In this summation, subindices should be taken modulo L, i.e., k L+ 1 should 
be understood as k\ and k L + 2 should be interpreted as k<i . 

Finally, consider the requirement that V(t) should be decreasing with 
time. For fixed discrete state i(t) = k and associated parameters a(t) = aj k 
and /3(t) = f3j k , the time derivative of the Lyapunov functional is 

^ V(t ) = 2 x(t) T P(A k x(t) + a k )-\- (8.13) 

+ 2 (a jk f T x(t) + f3 jk )f T (A k x(t) + a k ) (8.14) 



Here, j denotes the previous discrete state. The time-derivative is quadratic 
in x(t), and we can use the S-procedure to derive conditions for its negativity 
for all x that satisfy the mode invariant (similar to what was done in Chap- 
ter 4). This condition has to be imposed for all admissible pairs (j : k). We 
summarize the developments in the following result, which is a simplified 
version of the results in [73, 71] adopted to the notation of this book. 



Proposition 8.2— Lagrange Stability of Simple LHDS 

Consider the linear hybrid dynamical system (8.1) with parallel transition 

surfaces (8.9). For each admissible discrete transition j — > k, define 
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If there exists P = P T , Wjk >l 0, and (3jk such that 

P > 0 

A'kPjk + PjkAk + EfrWjkEk < 0 V(j, k) G T 

and o-jk, (3 jk satisfy (8.12) and (8.11) then the system is Lagrange stable. □ 

Note that the above formulation does not single out modes that contain 
equilibrium points and that the Lyapunov functional have affine terms in all 
modes. In order to verify the analysis conditions one thus needs to eliminate 
inherent equality constraints using the technique described in Section 4.12 
and solve the resulting optimization problem using a solver that supports 
equality constraints. 

The requirement that the integral term should be non-negative around 
every cycle of the transition diagram can first appear overwhelming since 
even very simple transition diagrams can admit a huge number of cycles. 
However, a fundamental result from graph theory (see, e.g ., [68, Chapter 4]) 
states that every cycle can be expressed as the sum of disjoint elementary 
cycles 1 . Thus, we only need to impose (8.12) for elementary cycles. 

The following example, taken from [71], illustrates the approach on a 
simple system that cannot be analyzed using Theorem 8.1. 

Example 8.3— Lagrange Stability of Hysteresis System 
Consider a linear system in a feedback connection with a hysteresis element 
as shown in Figure 8.5. Note that the hysteresis element is non-standard, in 
the sense that the hysteresis loop is clock-wise. We use the numerical values 
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Figure 8.5 Linear system in feedback interconnection with hystersis element. 

1 An elementary cycle is a cycle that does not visit any node (except for the initial and ter- 
minal nodes) more than once. The elementary cycles of a directed graph can be computed in 
polynomial time using, for example, the algorithm described in [156]. 
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The transition diagram for the linear hybrid dynamical system is shown in 
Figure 8.6, where we have introduced 



A\ = A4 = A 
Ei = 



A 2 = As = 



A B 







' c 


-1 






-c 1 


e 2 = e 4 = 


-c 
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E3 = 


C -2 



Since \f T x\ can become unbounded in modes 1 and 3, we require that 



OL 21 ^ 0 OL 41 > 0 U 23 > 0 0^43 > 0 



The transition diagram shows three fundamental cycles: 

£1 = (1,2,1) £2 = (3,4, 3) £3 = (1,2, 3, 4,1) 

However, since gu = 921 and # 34 = # 43 , the net amount of integral term is 
zero around the first two cycles. For the third cycle, noting that g 2 s = £34 
and # 4 i = 012 , the condition that the net amount of integral around the cycle 
should be non-negative can be simplified to 



1.6(0:34 - <*12) + (#3 4 - P12) > 0 



The LMI conditions of Theorem 8.1 have a feasible solution 



1.1663 

-1.4149 



-1.4149 

13.7400 



and 



ex 21 = 0:41 = ol 23 = a 43 = 0.0359 0:12 = —0.4922 U34 = —0.1925 

foi = Pai = 11.4197 = ^43 = 0 (3 12 = 0.0152 /?34 = 12.6413 

which guarantees Lagrange stability of the system. □ 



8.4 Comments and References 

Analysis of Linear Hybrid Dynamical Systems 

The area of hybrid systems has attracted a large interest over the last couple 
of years, and many tools and techniques for analysis and design of hybrid 
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xi = 2 




Figure 8.6 Transition diagram for linear system with hysteresis component. 



systems have appeared alongside with strong theoretical developments. A 
complete survey of the field is beyond the scope of this book, but relatively 
broad expositions can be found in, for example, [201, 135, 158, 52, 4, 5, 121] 

We have focused on stability analysis of a particular class of hybrid sys- 
tems, which we call linear hybrid dynamical systems. In the computer sci- 
ence literature, a linear hybrid system usually refers to a hybrid dynamical 
system in which trajectories of the continuous state in every mode are lin- 
ear (the continuous dynamics is described by integrators). Our definition is 
much more general. 

Linear hybrid dynamical systems have been studied by several research 
groups, see e.g., [100, 157, 76, 92]. Our main contribution has been to show 
how convex optimization can be used to construct piecewise quadratic Lya- 
punov functions for such systems. This allows computationally efficient anal- 
ysis of an important class of hybrid systems. The results presented in this 
chapter first appeared in [100, 102] and similar ideas have been suggested 
in the independent work [157]. The Lyapunov functional approach was in- 
troduced in [73, 71]. 

LHDS as Hybrid Automota 

A general framework for representing autonomous hybrid systems that has 
gained widespread acceptance is the hybrid automoton, see, e.g., [201]. Many 
variations of the basic model exist, but the following definition is useful to 
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enable comparisons between linear hybrid dynamical systems and other hy- 
brid system models. A hybrid automaton H is a collection 

H = (X, J, /, inv, jump, init) 



where 

X is the finite ordered set X = {aq, . . . , x n } of continuous state variables. 

I is the finite set of discrete modes. 

/ is a function / : R n x / i-> R n that defines the continuous flow in each 
discrete mode i e I through a differential equation x = /(x, i). 

inv is a set, inv c R n x /, called the invariant condition that describes the 
possible states ( z , i) of the hybrid automaton. For each discrete mode 
i, the set in v(i) = {z | (z, i) e inv} is called the mode invariant. 

jump is a set-valued function jump : R n x / i— > P(R n x /) called the jump 
condition . Here P(fi) denotes the power set of i.e., the set of all sub- 
sets of Q. The jump condition specifies if a transition from one discrete 
mode to another is possible, and what new value sould be assigned to 
the continuous variable after the jump. It is convenient to separate the 
jump map into guard conditions 

G (i(t ~ ) , i(t)) = {z(t~) G M n | 3z(t) G R n , (i(t), z(t)) G jump(^(i _ ), i(t - ))} 
and associated jump maps 

i(t), z(t~)) = {z(t) GM n | (z(t),i(t)) G jump(^(/: _ ), i(t - ))} 

init is a set init c R n x /, called the initial condition, that defines the possible 
initial states of H. 

We can now discern several features of linear hybrid dynamical systems. 
First, in each mode i the vector fields are affine in the continuous state vec- 
tor x. Second, the mode invariants are polyhedral regions in the continu- 
ous state space. Third, the jump maps do not alter the continuous state, i.e., 
i(t ), z(t~)) = z, and the guards are hyperplanes in the continuous 
state space. 
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9 

Concluding Remarks 



This book has treated analysis and design of piecewise linear control sys- 
tems. The developments contain two important ingredients: a powerful sys- 
tems class and a set of computationally efficient tools for system analysis 
and controller design. 

Piecewise linear systems is a natural extension of linear systems and a 
powerful model class for approximating smooth nonlinear systems. Many 
of the most common nonlinearities in engineering systems, such as satura- 
tions and relays, are piecewise linear. Although a large body of literature ex- 
ists on analysis of linear systems interconnected with specific components, 
almost no results exist for general piecewise linear systems. Research in the 
circuit and systems community has been focused on efficient simulation and 
static analysis, but have left analysis of the dynamics largely unattended. We 
have tried to take a broad view of piecewise linear dynamical systems and 
derived computational analysis tools that are generally applicable. 

Our description of piecewise linear systems draws from many sources 
but it also contains many new developments. We presented a simple matrix 
representation of polyhedral piecewise linear systems that is convenient for 
computations. We defined solution concepts and discussed issues related 
to attractive sliding modes. We extended uncertainty models for linear sys- 
tems to piecewise linear systems and demonstrated how this allows for less 
conservative analysis of uncertain nonlinear systems. We also showed how 
series, parallel, and feedback interconnections of piecewise linear systems 
are themselves piecewise linear systems. Such properties allow us to con- 
struct complex piecewise linear systems from simpler components and open 
up many interesting possibilities for input-output analysis of piecewise lin- 
ear systems. Finally, we reviewed memory-efficient representations of piece- 
wise linear systems with continuous right-hand side. 

We have provided a wide variety of analysis and design tools that are 
applicable to general piecewise linear systems. This includes simple proce- 
dures for computing equilibrium points equilibrium points and verifying in- 
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variant sets, but also sophisticated methods for stability analysis, gain com- 
putations, and controller design. The main contribution of the book has been 
to show how piecewise quadratic Lyapunov functions can be computed 
via convex optimization. This is a novel idea that can be implemented us- 
ing very efficient computations. The approach performs strictly better than 
analysis procedures based on the Popov criterion but its main attraction is 
that it easily deals with multi-variable nonlinearities. Moreover, the piece- 
wise quadratic Lyapunov functions can also be used in dissipativity analysis 
and the design of optimal control laws. An important feature of the analysis 
and design procedures in this book is that they provide simultaneous up- 
per and lower bounds on the estimated performance. These bounds can be 
used to judge the quality of the performance estimate and can also be use- 
ful in guiding partition refinements that aim at increasing the flexibility of 
the Lyapunov function candidate. When the piecewise linear system comes 
from the approximation of an asymptotically stable smooth nonlinear sys- 
tem, we have shown that such partition refinements will eventually allow 
us to compute a Lyapunov function using the suggested techniques. 

The basic analysis techniques have been extended in many directions. 
We have shown how smooth nonlinear systems and fuzzy logic systems can 
be analyzed rigorously in the piecewise linear framework and derived spe- 
cialized procedures for estimating regions of attractions. We also considered 
piecewise linear Lyapunov functions, showed how these can be computed 
via linear programming and suggested a procedure for automatic partition 
refinements based on linear programming duality. Finally, we extended the 
stability analysis procedures to a class of hybrid systems with piecewise lin- 
ear dynamics. These systems, which we have called linear hybrid dynami- 
cal systems, can be viewed as piecewise linear systems whose cells overlap 
on R n . Changes in the cell index are no longer determined directly by the 
continuous state, but have to be specified via a set of transition rules that 
depend on both the current cell index and the continuous state vector. 

Although this book has made significant progress towards a useful the- 
ory for piecewise linear systems, many problems remain open and many 
techniques can be improved. One of the strengths of this book is that it does 
not constrain the class of piecewise linear systems (e.g., to have continuous 
right-hand side) but presents procedures that are generally applicable. At 
the same time, this generality is also one of its weaknesses: it is hard to ana- 
lyze well-posedness, controllability, observability, etc. for general piecewise 
linear systems. The results presented here would benefit from a more com- 
plete system-theoretic understanding, even in a simplified and restricted 
setting. In addition, we believe that there is much to do on the computational 
side. By exploiting problem structure and developing specialized optimiza- 
tion methods, we believe that it would be possible to leverage the tools to 
allow analysis of complex systems with large state-spaces and partitions. 



143 




A 

Computational Issues 



The results derived in the previous chapters allow analysis and design of 
piecewise linear control systems based on numerical computations. For clar- 
ity of presentation, detailed discussions about computational issues have 
been postponed to this chapter. 

While most readers of this book probably have a good knowledge about 
linear and nonlinear systems theory [110, 176, 189, 116] they may be less fa- 
miliar with convex sets and convex optimization. Since these concepts play 
a central role in our approach, we will begin this chapter by an elementary 
overview of the material on convex sets and convex optimization that is 
used in this book. The presentation is inspired by the excellent survey [92, 
Chapter 1 1] but has been both restricted and extended to suit our purposes. 

Our computational analysis procedures have made heavy use of the con- 
straint matrices G it Ei and Fi without actually going into details of how 
these can be determined for a given partition. Using the results presented in 
the first part of this chapter we will show how these matrices can be com- 
puted for two important classes of partitions. In both cases, the constraint 
matrices can be computed efficiently using simple manipulations. We will 
also demonstrate how constraint matrices for complex partitions can some- 
times be computed from the descriptions of simpler partitions. 

The 5-procedure has played an important role in many computations 
and the conservatism of our approach is largely dependent on the conser- 
vatism of the 5-procedure. Our final effort will be to provide some further 
insight into the role of the 5-procedure and to derive two interesting results. 
Both consider simplex partitions. The first result states that analysis using 
the polyhedral 5-procedure is always less conservative than the use of min- 
imum volume ellipsoids. The second result establishes non-conservatism of 
the polyhedral 5-procedure for simplex cells in R n with n < 3. 
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A.1 Linear and Semidefinite Programming 

Convex optimization is the problem of minimizing a convex function over a 
convex set. Such problems have many desirable properties ( e.g ., locally op- 
timal solutions are also globally optimal) and admit an extensive and use- 
ful theory Formulating design problems as convex optimization problems 
is also practically appealing since these problems can be solved numeri- 
cally with great efficiency. For many classes of convex optimization prob- 
lems such as linear programming, semidefinite-programming, etc., public 
and commercial solvers exist that can exploit problem structure and solve 
large problem instances very efficiently. 

Very little of the general theory for convex optimization is needed in or- 
der to understand and apply the results of this book. All problems can be 
formulated as linear or semidefinite programming problems using classical 
linear algebra techniques and a handful of key results that are summarized 
below. Excellent presentations of the general theory and algorithms for con- 
vex optimization can be found in, e.g., [204, 210, 19, 28]. 

Linear Programming 

A linear program is a special class of convex optimization problems where 
the objective function and the constraints are linear 

minimize c T x 

subject to afx <bi i m 1 , . . . , m. 

The field of linear programming is well-developed, and large-scale lin- 
ear programs can be solved very efficiently using a variety of public and 
commercial codes (e.g., [46, 177, 6]). The theoretical aspects of linear pro- 
gramming are also well understood (see, e.g., [182, 21]). The only result that 
we will need, however, is Farkas’ lemma, which allows us to verify that a 
set of linear inequalities does not admit any solution. 

Proposition A. 1— Farkas’ Lemma 
Either there exists a vector x such that 

Ax A b 

or there exists a vector y hO such that 

y T A = 0 y T b < 0 



but not both. 



□ 
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Semidefinite Programming 

A semidefinite program is a convex optimization problem with a linear ob- 
jective and constraints specified in terms of linear matrix inequalities. A lin- 
ear matrix inequality (LMI) is an inequality on the form 

m 

F(x) = F 0 + J2 x iFi >0 (A.l) 

i= 1 

where x = (aq, . . . , x m ) is the variable, Fi = Fj e R nxn are given symmet- 
ric matrices and the matrix inequality F(x) >0 means that F(x) is positive 
semidefinite. Hence, a semidefinite programming (SDP) problem takes the form 



minimize c T x 
subject to F(x) > 0 



In this book, we will never write LMIs on the standard form (A. 1). Rather, 
we will phrase LMI constraints directly using matrix variables. For example, 
the constraint that the symmetric matrix variable X = X T e R 2x2 be posi- 
tive semidefinite can be written in standard form via 




Xu 


®12 


®12 


^22 




0 

0 



+ ^12 



0 

1 



1 

0 



+ ^22 



0 

0 



0 

1 



> o 



Clearly, any expression on the form 



Ao + 

3 i 



in an LMI in the symmetric matrix variables Xj and can be put in standard 
form by first expanding the matrix variables as above, and then pre- and 
post-multiply the basis matrices with the data matrices and Bij. Today, 
most semidefinite programming solvers (e.g , [62, 29, 122]) have good inter- 
faces that allow the user to specify LMIs directly in terms of matrix variables. 

A useful extension of SDP problems are the max-det problems 



minimize c T x + log det G (x) 1 
subject to G(x) > 0 
F(x) > 0 

which we will encounter when computing optimal ellipsoidal approxima- 
tions of polytopes. As of today, there is only a few codes (e.g, [206]) that can 
solve max-det problems. 

The following results are useful when formulating and manipulating lin- 
ear matrix inequalities (see, e.g., [27]). 
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Proposition A.2— Schur Complement 

Let Q(x) = Q(x) T , R(x) = R(x) T and S(x) be matrices whose entries de- 
pend affinely on the vector x. Then, the LMI 



Q(x) S(x) 
S(x) T R(x) 



is equivalent to 

R(x) > 0, Q{x) — S(x)R(x)~ 1 S(x) > 0 



□ 



Proposition A.3— S-Procedure [216, 27] 

Let Qo, . . . , Q p be quadratic functions of x e M n , 

Qi(x) = x T AiX + 2 bjx + Ci, i = 0, . . . : p (A.2) 

where A* = Af. Consider the following condition 

Qo(x) > 0 for all x such that Qi(x) >0, i = 1, . . . ,p (A.3) 

A sufficient condition for (A.3) to hold is that there exist positive scalars 
Ti > 0 such that 



p 

Qo(x) > J2n Ql (x) (A.4) 

i= 1 

Moreover, if p = 1 and there exists an xo such that Qo(^o) > 0, then condi- 
tion (A.4) is also necessary. □ 

When using the S-procedure in this book, we will re-write the quadratic 
inequalities (A.2) in matrix format 



X 


T 


hi? 




X 


1 




o 

hO 




1 



and then convert (A.4) into an LMI condition using the following result. 
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Proposition A. 4 

Let Ai = Af € R nxn , bi e R n and c* e R. Then, 



if and only if 



X 


T 


-Ai bi 




X 


1 








1 



> 0 Vx e R n 



bf 



bi 



Ci 



> 0 



(A. 5) 



(A. 6) 

□ 



Note that the when we replace the non-strict inequalities in Proposition A.3 
and Proposition A.4 by strict inequalities, Condition (A.4) is only sufficient 
for (A.3) and Condition (A.6) is only sufficient for for (A.5). For example, the 
inequality 



X 


T 


0 


o" 




X 


1 




0 


1 




1 



VxeR n 



cannot be verified using an analogue of Proposition A.4 that uses strict in- 
equalities, since the matrix 



0 0 
0 1 

is only positive semidefinite. However, in our experience, the conversion 
from a strict quadratic inequality to a strict LMI and the use of the strict S- 
procedure appears to work very well in practice. A more detailed discussion 
about the specific way in which the 5-procedure is used in this book will be 
given in Section A.5. 



A.2 Polyhedra, Polytopes, and Ellipsoids 

As stressed in Chapter 2, any analysis procedure for piecewise linear sys- 
tems must adequately account for both the partition and the dynamics. In 
this section, we will review the necessary background material on describ- 
ing and manipulating polyhedra, polytopes and ellipsoids. 

The notation in this section is standard with (perhaps) a few exceptions: 
for a matrix M, Mj_ denotes any full rank matrix such that Mj_M = 0; Mjsf 
is a full rank matrix whose columns form a basis for the null space of M\ 
and Aft denotes the generalized (or Moore-Penrose) inverse of Af . 



148 




A.2 Polyhedra, Polytopes , and Ellipsoids 



Convex Sets and Their Representations 

Convex sets can be represented in many ways. In this book, we make use of 
two different representations: the constraint representation 

C = {x G R n | g(x) < 0 A h{x) = 0} 

and the constrained parameter representation 

C = {x = f{9) | g(6) < 0 A h{9) = 0 9 G R m } 
where g(-) is convex and h(-) is linear. 

Hyperplanes and Halfspaces 

Definition A. 1— Hyperplane 
A hyperplane in R n is a set on the form 

dH = {x G M n | a T (x — xo) = 0} (A. 7) 

It can be equivalently represented as 

dH={x = x 0 PAO | OeR 71 - 1 } (A.8) 

□ 

Geometrically, the parameter a is the normal vector of the hyperplane, and 
xq is an arbitrary point on the hyperplane, see Figure A. 1 (left). 



on 





Figure A.l Hyperplane (left) and halfspace (right). 
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To convert from the constraint representation (A. 7) to the constrained 
parameter representation (A.8) we let A = ajy, to pass from (A.8) to (A. 7) 
we let a T = A±. In some cases, the hyperplanes will be given on the form 

<97 i = {x G R n | a T x = b} 

We can convert from this representation to (A. 7) by setting xo = (a T )f 6. 

Definition A. 2— Halfspace 
A halfspace is a set on the form 



H = {x g R n | a T (x — x o) < 0} 



(A. 9) 



It can be equivalently represented as 

H = {x = xo + AO + at | 0 g R n_1 , t > 0} (A. 10) 

□ 

To pass from (A. 9) to (A. 10) we let A = ajj. To convert from the representa- 
tion (A. 10) to (A.9) we let a T = -sign(A^a)A_L. 

Polyhedra and Polytopes 

Definition A. 3— Polyhedron 
A polyhedron is a set on the form 

V = {x G R n | ajx < bi, i = l,...,m} = {xG M n | Ax ^ b} (A.l 1) 
It can be equivalently represented as 

V = {x = V\yW!i \ = A^O, /x^O} (A. 12) 

i 

□ 

The representation (A. 11) describes the polyhedron as the intersection of 
a finite number of halfspaces while the formulation (A. 12) represents the 
polyhedron as the set addition (Minkowski sum) of the convex hull of the 
columns of the matrix V and the conic hull of the columns of the matrix 
W. The columns of V are called the vertices of the polyhedron, while the 
columns of W are called extreme rays . 
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Definition A. 4— Polytope 

A polytope is a bounded polyhedron. It can be represented as 

V={x\ Ax<b} (A. 13) 



or, equivalently, 



p = { x = yA| ^A i = l, XhO} (A. 14) 

i 

□ 

The representation (A. 13) is sometimes referred to as the halfspace representa- 
tion, while (A. 14) is called the vertex representation. The conversion between 
the two representations (A. 13) and (A. 14) is conceptually simple but com- 
putationally hard (see, e.g., [9]). If scalability and generality of an algorithm 
is important it is usually advised to use only one of the representations and 
not introduce steps that require the conversion between the two. For small- 
to moderate-sized problems, the conversion can be made using a variety of 
publically available codes, see [41,9] and the references therein. 

The dimension of a polyhedron is the dimension of its affine hull. For ex- 
ample, the dimension of a singleton {v} is 0, the dimension of a line segment 
is 1, etc. 

Definition A. 5— Faces, Facets, Vertices and Edges 
L et V be a polyhedron and let H be a halfspace such that V C H and V F = 
V fl H ^ 0. The polyhedron V F is called a face of V. Faces of dimension 0 
are called vertices, faces of dimension 1 are called edges, and faces of V of 
dimension dim(V) — 1 are called facets, see Figure A.2. □ 

We note that polytopes defined by few halfspaces might have very many 
vertices and vice versa (see, e.g., [9]). As a simple example, note that a unit 
cube in R n is defined by 2 n halfspaces but has 2 n vertices. 

Two particular families of polytopes, simplices and parallelotopes, will 
appear repeatedly in this book. They are defined as follows. 

Example A. 1— Simplex 

An r-simplex in R n is a polytope that can be expressed as 

( r+ 1 

V = X = ^2 A iVi, ^2 Xi = 1, A i 

K. i= 1 i 

for affinely independent vectors vi,...,v r . □ 
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vertex 

Figure A. 2 Facets, edges, and vertices of a polytope. 



Example A. 2— Pakallelotope 

An r-parallelotope in R" is a polyhedron that can be expressed as 
V {.r e M” | \4'x - t H \ < 1 , * 1 , . . . , r} 

for linearly independent vectors a \ . . . . , a,,.. □ 



Quadratic Sets and Ellipsoids 

A natural extension of polyhedra, which are defined by linear inequalities, 
is to consider sets defined by quadratic inequalities in the following way 

Definition A. 6— Quadratic Set 
A quadratic set is a set on the form 

Q = {x e R n I x r Ax + 2b 1 ' x + c < 0} (A.15) 

□ 



It is often convenient to parameterize the quadratic polynomials as 



X 


T 


'a b 




X 


1 




b T c 




1 



where .4 A 1 ' g R” xr \ b G R n and c € R* 

Quadratic sets can be convex or non-con vex, bounded or unbounded, 
connected or disconnected (see, e.g [92])* Of particular interest are ellip- 
soids since they define bounded, connected and convex sets. 
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Definition A. 7— Ellipsoid 
A n ellipsoid is a set on the form 



A.2 Polyhedra, Polytopes , and Ellipsoids 



S = { x \{x — x 0 ) T P 1 (x — x 0 ) < l} (A. 16) 

with P = P T > 0. It can be equivalently represented as 

£ = {x = Q0yx 0 | ||0|| 2 < 1} (A. 17) 

□ 

Geometrically, xo is the center of the ellipsoid while the matrix P determines 
its shape. To be more specific, let P = UAU T be an orthonormal eigenvalue 
decomposition of P, where U = [u\ ... u n \ and A = diag(Ai, . . . , A n ) with 
Ai > 0 (such a decomposition is always possible when P is positive definite). 
Then, the semi-axes of £ are given by ui and their respective lengths are y/Al, 
see Figure A. 3. 




Figure A.3 Center, semi-axes, and semi-axes lengths for an ellipsoid. 

We can convert from the representation (A. 16) to (A. 17) by letting Q = 
P 1 / 2 . If the ellipsoid is defined by a set on the form (4.30), then we can pass 
to the representation (A. 16) by letting P = (b T A~ x b - c)A~ 1 f x$ = -A -1 b. 
We will also encounter ellipsoids parameterized as 

£ = {xeR n I \\Cxyd\\ 2 < 1} 

We can pass from the representation (A. 17) to this representation by letting 
C = Q- 1 and d = -Q _1 x 0 . 
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In many cases, it is desirable to optimize the “size” of an ellipsoid subject 
to constraints. A natural measure of size is the volume of the ellipsoid, 

vol(5) = v n V det P 

where v n is the volume of the unit sphere in R n . Other measures of size in- 
cludes the trace of P (which measures the sum of squared semi-axes lengths) 
and the maximum eigenvalue of P (which measures the maximum semi- 
axis length). 

Containment 

A fundamental problem that arises on several occasions in this book (e.g., 
Section 7.1) is to determine if a convex set contains another convex set. The 
following two results are then useful. 

Proposition A. 5— Ellipsoid in Halfspace 

The ellipsoid (A. 16) is contained in the halfspace (A. 9) if and only if 

a T x o + V a T Pa < b (A. 18) 

If the ellipsoid is represented as (A. 17), the containment condition is 

o T x 0 + ||Qa|| 2 < b 



□ 

To verify that an ellipsoid 8 is contained in a polyhedron V, the containment 
condition (A. 18) must be satisfied for all halfspaces that define V. 

Proposition A. 6— Quadratic Set in Quadratic Set 
L et 



Qi = {x G R n | x T AiX + 2 bjx + q < 0 } i — 1, 2. 

be two quadratic sets and assume that Q 2 ^ 0- Then Q\ c 0,2 if and only if 

A 2 b<2 
}>2 c 2 

for some positive scalar r > 0. □ 



At h 
bi ci 



< 0 



154 
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Ellipsoidal Approximation of Polytopes 

Although we consider polyhedral piecewise linear systems, we have seen 
that it can sometimes can be useful to use (approximate) ellipsoidal descrip- 
tions of the polyhedral cells. For example, in Section 4.9 we demonstrated 
how ellipsoidal cell boundings allow for faster computations of piecewise 
quadratic Lyapunov functions, and in Section 6.1 we showed how the use 
of ellipsoidal cell descriptions enables a convex formulation of the quadratic 
stabilizability problem for piecewise linear systems. The following results 
demonstrate how optimal ellipsoidal approximations of polytopes can be 
computed using convex optimization. 



Proposition A. 7— Ellipsoidal Approximation of a Polytope 
Consider the polytope 



V=co 

The parameters (C, d) that specify the minimum volume ellipsoid 
fmve = {^ G M n | || Cx + C?|| 2 ^ 1} 
that contains V can be obtained by solving the max-det problem 

minimize log det C ~ 1 

I Cvi + d] 

vJC T + d T 1 

c = c T > 0 



subject to 



> 0 



= 1 , 



Moreover, the ellipsoid 

£' = {xeR n I \\n(Cxyd)\\ 2 <1} 

is contained in V. □ 

The inner ellipsoid obtained by scaling the minimum volume ellipsoid is of- 
ten called the Lowner-John ellipsoid and the scale factor n can be improved 
to y/n if the polytope is symmetric around its center [105]. When V is a sim- 
plex, we can obtain the minimum volume ellipsoid directly. 

Proposition A. 8— Ellipsoidal Approximation of a Simplex 

Let X be an n-simplex in R n with vertices v it and let V = [v\, . . . ,v n +i]- 

Then, the minimum volume ellipsoid that contains X is given by 

S mve ={xeR n | X T (VV T )- 1 X < 1 } 

□ 

The analogue result of Proposition A. 7 for polytopes in halfspace represen- 
tation is the following. 
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Proposition A. 9— Ellipsoidal Approximation of a Polytope II 
Consider the polytope 

V = {x e R n \ ajx < bi,i = 1, , m} 



The parameters (Q, x$) that specify the maximum volume ellipsoid 
£mve = {x = Qu + xo | IMI2 < 1} 
contained in V are obtained by solving the max-det problem 



minimize log det Q 1 

\{h - afx 0 )I 



subject to 



( Qai) T 

Q = Q T > 0 



Qa, { 

a- x 0 



bj — rtT 



>0 



i — 1, . . . , m. 



Moreover, the ellipsoid 

S f = {x = nQu + xq | \\u\\2 < 1} 



contains V. □ 

Similar to above, the factor n can be replaced by yjn if the poly tope is sym- 
metric around its center. A drawback with this formulation is that the opti- 
mized ellipsoid is contained in V and that the scaled ellipsoid that contains 
V is not necessarily very tight. 

For the special case of parallelotopes, we can obtain an explicit expres- 
sion for the minimum volume ellipsoid that contains V. 



Proposition A. 10— Ellipsoidal Approximation of Parallelotope 
Let X be an n-parallelotope in R n , 



X = {x G R n | | af x — bi | < 1, i = 1, . . . , n} 



and define 



H = 



Then, the minimum volume ellipsoid that contains X is 

£mve = {x e R n | X T H T Hx < Ti } 



□ 
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A.3 Polyhedral Partitions 

More complex geometrical structures can be constructed by “gluing” to- 
gether simple polyhedra, see Figure A.4. In computational geometry and 
combinatorics, the resulting objects are called polyhedral complexes (see, e.g., 
[221]). In our setting it is more natural to adopt a dual view, since the com- 
plexes arise from the subdivision of a given set into convex polyhedra. We 
have called these objects polyhedral partitions, and defined them as follows: 

Definition A. 8— Polyhedral Partition 

A polyhedral partition X C R n is a collection of closed n-dimensional convex 
polyhedra Xi such that 

1. U ieI Xi = X 

2. the polyhedra have disjoint interior, int(X^) n int (Xj) = 0 if i ^ j 

3. the polyhedra only share their common boundaries. In other words, if 
Xi H Xj ^ 0, then Xi n Xj is a common face of Xi and Xj. 

□ 

In other fields of science, similar objects have been defined under the names 
piecewise linear manifolds [2] and, if all polyhedra have the same shape, meshes 
(see, e.g., [18]). 

In order to represent a polyhedral partition we need to represent the in- 
dividual polyhedra as well as their topology. The development of efficient 
data structures for polyhedral partitions is far from trivial and the complex- 
ity of the data structure depends on what information we need to be able 
to access quickly. For example, in some situations we only need to represent 
the n-dimensional polyhedra and their adjacency relations while in other 
cases we might need to access definitions and adjacency information of all 
cell faces (down to edges and vertices). Data structures for efficient represen- 
tation of polyhedral partitions can be found in, for example, [159, 36, 130]. 

As we have seen in Section A. 2, convex polytopes admit two equivalent 
representations: they can either be represented as the convex hull of a finite 
number of points or as the intersection of a finite number of halfspaces. From 
these two descriptions, two classes of partitions appear particularly natural: 
one class is the partitions induced by a number of hyperplanes and the other 
is the partitions induced by a number of points. 



Hyperplane Partitions 

A hyperplane partition is a partition which is induced by a number of global 
hyperplanes. The cells of a hyperplane partition are convex polyhedra that 
have these hyperplanes as boundaries. Given a set of convex polyhedra, 
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an associated hyperplane partition can be obtained by extending the cell 
boundaries globally. This extension of boundaries may induce new cells, 
see Figure A.4. Hence, a partition obtained in this way may have more cells 
than the number of polyhedra that generated it. 





Figure A.4 A hyperplane partition generated from an initial set of polyhedra by ex 
lending Llieir boundaries globally. 



Simplex Partitions 

A simplex partition is a partition induced by a number of grid points. Given 
a set of grid points in R n , we can form a simplex partition by creating n- 
simplices that contain exactly n + I grid points and form a partition. Note 
that a set of points in R n can be combined into many diff erent simplex par- 
titions, see Figure A. 5. In other words, a set of vertices does not induce a 
unique simplex partition (see 1 124| for further details). 

Every convex polytope can be partitioned into simplices by the possible 
insertion of new vertices. Thus, given some initial polytopic partition, it can 
always be refined into a simplex partition. A simplex partition derived in 
this way may then have more cells than the original (non-simplex) partition 
that generated it, see Figure A, 5. 



V X l/ 2 V 1 




Figure A. 5 Two simplex partitions induced by the four vertices i'k of a square. 
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A.4 Constraint Matrices 



The analysis procedures developed in this book are not restricted to any par- 
ticular partition type but their formulation do require that some key prop- 
erties of the partition are expressed in the specific matrix format defined in 
Section 4.5. The matrices used in different computations are the cell iden- 
tifiers, G it the cell boundings, E it and the continuity matrices, Fi . For this 
parameterization to be useful we need to be able to compute these matrices 
for a given partition. In this section, we will show how this can be done for 
hyperplane and simplex partitions. 

Cell Identifiers and Cell Boundings 

The cell identifiers Gi and cell boundings Ei are closely related. In fact, the 
only reason that the cell identifiers cannot be used in the Lyapunov func- 
tion computations is that they typically do not satisfy the zero interpolation 
property which is critical for formulating the Lyapunov function search us- 
ing strict LMIs. The following procedure shows how cell boundings with 
the zero interpolation property can be computed from the corresponding 
cell identifiers. 



Algorithm A. 1— From Cell Identifier to Cell Bounding 

Let {Xi} ieI be a polyhedral partition with associated cell identifiers Gi. The 

corresponding cell boundings can be computed as follows: 

if i e Io, then Ei is obtained by deleting all rows of Gi whose last entry 
is non-zero. 



if i e I\ and Xi is unbounded, then Ei is obtained by augmenting Gi 
with the row [0i xn 1], otherwise Ei = Gi. 



□ 



It is natural to ask if there is any loss in using cell boundings rather than 
cell identifiers in conditional analysis of piecewise quadratic functions. The 
following lemma shows that this is not the case. 



Lemma A.l 

Consider the piecewise quadratic function V{x) constructed in Lemma 4.2 
and let Ei be cell boundings constructed using Algorithm A.L Then Ei have 
the zero interpolation property, and 

V(x)>0 fovxeXi (A. 19) 

if and only if 

V{x) > 0 for {x | EiX h 0} (A.20) 

□ 
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Proof: See Section B.6. 



Since the cell boundings can be computed from the cell identifiers with- 
out introducing any conservatism we only need to derive procedures for 
computing cell identifiers and continuity matrices. We will now show how 
this can be done for hyperplane and simplex partitions. 



Constraint Matrices for Hyperplane Partitions 

A hyperplane partition is a partition induced by K hyperplanes, 

dHk = {x | h%x + g k = 0} k=l,...,K. 

The partition induced by a saturated linear state feedback, illustrated in Fig- 
ure 2.1, is a typical example. For convenient representation, we collect all 
hyperplane data in a hyperplane matrix H : 



hj 



H = 



91 



L 9 k \ 

We adopt the convention that every hyperplane is defined with g k < 0. Each 
hyperplane induces two closed halfspaces, 



Hi = {x | hjx + g k > 0} 
Kk ={x | hlx + g k < 0} 



which we will call the positive and negative induced halfspace of dH k , re- 
spectively. The convention g k < 0 then implies that 0 e Ti j~ for all k. 

Every cell of the partition can be specified by stating whether it belongs 
to the positive or negative induced halfspace of each hyperplane. Hence, a 
cell identifier Gi for cell Xi can be obtained by multiplying the kth row of 
Ft with -1 if Xi C H ^ and with +1 if Xi C Note that the resulting cell 
identifier might have many redundant constraints. 

Continuity matrices F* can be computed as follows. Let the &th row of 
Fi be equal to the kth row of H if Xi C H £ and equal to the zero vector 
otherwise. Since 0 € H ~j~ for all k , this assures that 

Fix = max { Fix , 0} x G X iy 

where ma x(z, v) denotes element-wise maximum. This implies that the ma- 
trices Fi have the zero interpolation property. We summarize the develop- 
ment in the following proposition. 
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Proposition 

A. 1 1— Constraint Matrices for Hyperplane Partitions 
Let {Xi} ie i be a hyperplane partition. The matrices G* and Fi constructed 
as above satisfy the conditions (2.4) and (4.14), respectively. Moreover, the 
matrices Fi have the zero interpolation property. □ 

In order to give the continuity matrices full column rank, we can augment 
them according to 



p Ft fi 

F, = I 0 



(A. 21) 



The constraint matrices for the saturated system given in Example 2.3, Ex- 
ample 4.4 and Example 4.5 were computed using the procedure outlined 
above and augmented as in (A.21). The following example illustrates the 
use of the hyperplane matrix when computing the cell identifier for one of 
the regions. 

Example A. 3— Constraint Matrices for Saturated System 
Consider again the linear system with actuator saturation, 

x = Ax + bsat(k T x) 

The hyperplanes induced by the saturation give the hyperplane matrix 



hj gi 




' k T -l' 


h-2 92_ 




-k T -1 



Note that gk < 0 for all k. Consider the cell X 3 corresponding to positive 
saturation (k T x > 1). Since X 3 C f-i\ and X 3 C we can obtain a cell 
identifier by multiplying the second row of H with -1. This gives 



£3 



k T -1 
k T 1 



As only dHi is a boundary of X 3 , we delete the second row and arrive at 



G s = 



k T 




which was the cell identifier given in Example 2.3. 



□ 
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Constraint Matrices for Simplex Partitions 

A simplex partition X C R n is induced by a set of points {z/o, . . . , vk }• We 
will assume that the origin is one of the vertices and, without loss of gener- 
ality, let v o = 0. Each cell Xi is defined as the convex hull of n + 1 affinely 
independent points from this set in such a way that {Xi} ie i form a polyhe- 
dral partition of co{z/q, . . . , z/#}, see Figure A. 6. 




Figure A.6 Simplex partition of state space. 

Since cells are represented by the convex hull of their vertices, each x e 
Xi can be represented as 

x = z k v k x e Xi (A.22) 

k : j/feEX* 

with z k > 0 , Zk te 1 - The numbers z k are sometimes called the Barycentric 
coordinates of X{. Since each cell has n \ 1 affinely independent vertices, the 
Barycentric coordinates are unique and can be obtained by solving 

f X =Ys k VkZ k 
1 1 =£ k z k 

If the decomposition (A.22) is used for x ^ Xi then at least one of the 
Barycentric coordinates will be negative. This indicates that the affine map- 
ping from x to the Barycentric coordinates of Xi qualifies as cell identifier 
(this is indeed what we will use later). 

Similarly, each x e X can be uniquely represented by a convex combina- 
tion of all vertices of the partition 



K 

x = y^z k v k 
k = 0 



(A. 2 3) 
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by letting Zk = 0 if i/k & Xi. Clearly Zk > 0 for all k and = 1* Since 

the partition coordinates are unique for each x e X, the affine mapping 
from x to z = [zo ... zk] is continuous and can be used for constructing 
continuity matrices. 

To describe the computations we introduce an extraction matrix Si e 1 (tf+hxfa+i) 
as follows: the kth row of Si is zero for all k such that z/& Xi and the non- 
zero rows of Si are equal to the rows of an identity matrix. The extraction 
matrix has the property that 

z(x) = S{zi(x) x e Xi 

where Zi is the vector of Barycentric coordinates of x e Xi (the precise or- 
dering of the components of % depends on the choice of extraction matrix). 
Also introduce 



Then, we have 




(A. 24) 



x = Vz = VSiZi 

z = s^i = Si(ySi)- x x 

The matrix VSi is invertible since the vertices of each simplex cell are affinely 
independent. It turns out that continuity matrices constructed from the map- 
ping xi —> z does not have the zero interpolation property. However, such 
matrices can be constructed by simply disregarding the partition coordinate 
zo in the mapping. Thus, we propose to use 



Fi = 



0 I K 



Si{VSi 






which guarantees that 



FjX = 



0 I K 



Hi • •• z k 



(A. 25) 



Since the origin corresponds to zo = 1 and Zi = 0, i = 1, . . . , K in (A. 2 3), 
we have Fi [ 0 1] T = 0 so these matrices do indeed have the zero interpo- 
lation property. Cell boundings with the zero interpolation property can be 
constructed as 



Ei 



= Et 



0 

Fi 



We have the following result. 



(A. 26) 



163 




Appendix A. Computational Issues 

Proposition A. 12— Constraint Matrices for Simplex Partitions 
Let {Xi} ieI be a simplex partition. The matrices f* and Ei constructed as 
in (A. 25) and (A. 26) have the zero interpolation property and satisfy the 
conditions (4.14 and (4.19), respectively. □ 

It can be verified that the cell boundings computed above are equivalent to 
the ones that are obtained by first computing cell identifiers 

Gi = iVEi)- 1 

and then applying Algorithm A.l. 

The above construction extends straightforwardly to unbounded poly- 
hedra by considering simplices that have vertices “at infinity”. In this case 
every xGlj can be written as 



q K 

X = ^2z k lS k + ^2 Z k w k 

k = 0 k=q-\-l 

where z k > 0 and J2 k =o Zk = 1* The vectors i/ 1: . . . , v q are vertices of the 
polytope, while w q + 1 , . . . , w K define extreme rays. The computations of con- 
straint matrices remain the same and statements above hold true also in this 
case, provided that each cell has at least one vertex and that we define 



ISq . . . Vq kVq - (-1 • • • kUp 

1 ... 1 0 ... 0 



Building Complex Partitions from Simple Partitions 

In the computations above we assumed that the dimension of the partition 
was the same as the dimension of the state space, i.e., that the partitioning 
has been done with respect to all state variables. However, when significant 
nonlinearities are confined to a subset of the states (as in the min-max se- 
lector system described in Section 4.6) it can be natural to concentrate the 
flexibility of the Lyapunov function candidate to these states. 



Partitioning a Subset of the State Let x e M n and assume that the parti 
honing has been performed on the subspace 

Z = {z eR q \ z = Cx} 

Then, constraint matrices for the corresponding partition in R n can then be 
constructed by post-multiplying the matrices constructed in R q by 



N = 



C 

0 



0 

1 
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That is, let F Z i , E Z % and G Z i be constraint matrices for a polyhedral partition 
in R q , let C G R? x ( n+1 ), and let N be defined as above. Then, 



Fi = F Zi N : Ei = E Zi N , Gi = G Zi N , 

are constraint matrices for the corresponding cells in R n . Moreover, if F Z i 
and E Z i have the zero interpolation property, then so have Fi and This 
approach was used to construct constraint matrices for the min-max selector 
system analyzed in Chapter 4. If the constructed continuity matrix does not 
have full row rank, this can be achieved by the augmentation (A.21). 

Creating Cells by Intersecting Partitions Another issue appears when we 
interconnect two piecewise linear systems for which we have already com- 
puted constraint matrices. Let S x be a piecewise linear component with state 
vector x G R nx and partition {. Xi} ie i , and let S z be a piecewise linear com- 
ponent with state vector z G R n * and partition {Zj}j E j. The interconnected 
system can then be realized with a state vector v e R n *+ n *. The partition of 
the interconnected system is obtained as the product between the partitions 
of the components 



1 1 ^ i->j ^ jy 



where 



Fj,j — z'l | x G Xij z G Zj } 



The corresponding constraint matrices can be constructed by first extend- 
ing the constraint matrices for the subsystems into R nx x R n * (using the 
technique described above) and then stacking them on top of each other 
(creating the intersection). Similar to above, if the constraint matrices of the 
individual components have the zero interpolation property then so have 
the matrices describing the product partition. This approach was used to 
construct constraint matrices for the fuzzy system example in Section 7.4. 



A.5 On the S-Procedure in Piecewise Quadratic Analysis 

The 5-procedure plays a crucial role in the piecewise quadratic analysis, as 
it allows us to verify conditions on the form 

x T Pix > 0 for x G Xi 
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using semidefinite programming. The conservatism of our approach de- 
pends to a large extent on the conservatism of the 5-procedure. 

The 5-procedure has been studied in depth by Yakubovich in the early 
70’s [216]. To review his key results, let 

V (x) = x T Px 

be a quadratic function, and let 

Si = {x | x T Six >0} i = 1, . . . , m 

be quadratic sets. Yakubovich established that the condition 

V(x) > 0 Vx g Si 

holds if and only if there exists a non-negative scalar ui such that 

P — UiSi > 0 

Clearly, if there exist positive scalars Ui such that 



m 

P-^2uiSi> 0 (A. 2 7) 

i= 1 

then 



V(x) >0 Vx G P%Li£i (A. 28) 

By a simple example Yakubovich proved that the converse is, in general, not 
true. In other words, there exist quadratic functions V(x) and quadratic sets 
Si such that (A. 28) holds but where no solution ui > 0 to the LMI (A.27) can 
be found. Hence, while the 5-procedure is a necessary and sufficient condi- 
tion for positivity of a quadratic function on a quadratic set it is only suffi- 
cient for verifying positivity on the intersection of several quadratic sets. 

To understand the implications for our approach, let ef k denote the kth 
row of Ei. Then, our use of the 5-procedure is a special case of (A.27) since 

P - EfUiEi = P - Ujk&ije J k > 0 (A. 29) 

j,k 

This indicates that the 5-procedure is only a sufficient condition for verify- 
ing positivity of a quadratic function on a polyhedron. In this section, we 
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will show that our use of the 5-procedure is conservative in general, but 
that no conservatism is introduced for simplex partitions in R n with n < 3. 

It is easy to come to the premature conclusion that it is more conservative 
to use polyhedral cell boundings (for which our use of the S-procedure is 
only sufficient) than to fix a quadratic cell bounding of the polyhedral region 
and use this in the LMI computations (so that the S-procedure step itself 
becomes lossless). We will show that this is, in general, not true. For some 
important classes of partitions we will be able to prove that the polytopic 
relaxation is always stronger than quadratic relaxations. 

Copositivity and Non-conservatism of the S-Procedure 

The problem of verifying positivity of a quadratic form on a polyhedron 
has a long history with numerous applications in optimization (see, e.g., [48, 
167, 89]). Most of the work has been concerned with verification of matrix 
copositivity. A matrix C is said to be copositive if x T Cx > 0 for all x e R™ . 
For some time it was conjectured that if C is copositive, then it can be written 
as the sum of two matrices 



C = P y N (A. 30) 

where P is positive semidefinite and N has non-negative entries. The fol- 
lowing result was proven in [48] (see also [89]). 

Proposition A. 13— Decomposition of Copositive Matrices 

For dimensions n < 4, every copositive matrix C can be decomposed as in 

(A. 30). However, indecomposable copositive matrices exist for n > 5 . □ 

The decomposition (A. 30) suggests an immediate test for copositivity: if, for 
a given matrix C, there exists a matrix N with non-negative entries such that 

C-N> 0 (A. 31) 

then C is copositive. However, as stated in Proposition A. 13, this test is only 
sufficient and there exist copositive matrices for which the test fails. To see 
the relation to our approach, consider the problem of verifying positivity of 
x T PiX on the positive orthant, for which Qi = [I 0] and Ei = /. Then, the 
positivity condition (A. 29) reduces to 

Pi - Wi > 0 

which is identical to (A. 31). Hence, by Proposition A. 13, there exist quadratic 
functions that are positive on polyhedra but where this fact cannot be veri- 
fied using our approach. 
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Proposition A. 13 can also be used to prove that our approach is non- 
conservative in some cases. In what follows, let Xi be a simplex in R n with 
associated continuity matrix F* and cell bounding Ei computed as in Sec- 
tion A. 4. Assume that i e I\ and consider verification of the inequality 

x T FjTFiX >0 x £ Xi (A.32) 

Since i e Ii, we have Ei = Gi = {VECT 1 an d % = & \ x are the Barycentric 
coordinates of x e Xi. Condition (A.32)is equivalent to 



trsT 



0 I K 



0 I K 



EjZj > 0 



Zi G 



|( n +l) 



(A. 3 3) 



i.e. f verification of copositivity of % = Ej[0 I k ] t T[ 0 Ik]£% (the restriction 
1 T Zi = 1 can be disregarded due to homogeneity). Proposition A. 13 states 
that for n < 4, Ti is copositive if and only if there exists a Wi >z 0 such that 

Ti-Wi >0. (A. 34) 



This inequality holds if and only if 

(V£ i )- T T i (V£,;r 1 - (V£ i )- T W i (V£ i )~ 1 > 0 



which can be reformulated as 

PfTF - EfWiEi > 0. 

This is a non-strict version of the LMI condition in Theorem 4.1. The follow- 
ing proposition extends the result to strict inequalities and the exact formu- 
lation used in Theorem 4.1. 



Proposition A. 14— Non-conservatism of the ^-procedure 

Let {Xi} ie i be a simplex partition in R n with n < 3, and with constraint 

matrices Ei and Fi computed as in Section A.4. Then 

V(x) = x T F^TFiX > 0 for x G X^\{0}, i G I 



if and only if there exists matrices Wi with positive entries such that 

F?TFi - EfWXi > 0 for i e I 0 

FfTFi - EfWiEi > 0 for i e h 

□ 

Proof: See Section B.6. 
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Polyhedral Relaxation is Stronger than Ellipsoidal 

In Section 4.9, we illustrated by a simple example how the use of quadratic 
cell boundings can allow for significant computational savings compared to 
a formulation based on polyhedral cell boundings. However, as the same 
example demonstrated, these savings come at the price of increased conser- 
vatism in the analysis. The developments in Section A. 4 will now allow us to 
be more precise: we will show that if the piecewise quadratic computations 
that use the quadratic cell boundings of minimum volume have a solution 
then so do the computations based on polyhedral cell boundings, while the 
opposite is not true. This is contrary to a statement in [72, Section 4.1] where 
it was indicated that computations using ellipsoidal cell boundings would 
be less conservative than those using polyhedral relaxations since the S- 
procedure may be lossy when several quadratic terms are used. We have 
the following result. 

Proposition 

A. 15— Polyhedral Relaxation is Stronger than Ellipsoidal 
Let Xi be a simplex with associated cell bounding Ei (computed as de- 
scribed in Proposition A. 12 and Algorithm A.l), and let Si describe the min- 
imal volume ellipsoidal bounding computed using Proposition A.8. Then, if 
theLMI 



Pi-TiSi >0 (A. 35) 

has a solution then so has the LMI 

Pi - EfWiEi > 0, (A. 36) 

but there are cases when (A. 36) admits a solution while (A.35) does not. □ 
Proof: See Appendix B.6. 

To understand this result better, note that the role of the ^-procedure in 
the verification of the constraint 

V{x) = x T PiX >0 x G Xi (A. 37) 

is to compute a quadratic set that contains Xi but does not intersect the 
set Vf~ = {x | x T PiX < 0}. Clearly, the volume of the covering ellipsoid 
might have very little to do with this separation. Figure A.7 illustrates this 
for the specific counter-example used in the proof of Proposition A. 15. Here, 
the minimum volume ellipsoid that contains Xi intersects the set V~ and 
can therefore not be used to verify the desired inequality. The polyhedral 
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Figure A. 7 The counter example in Proposition A. 15. The minimal volume ellipsoid 
fails to separate Xi from V { ~ (left) , while optimizing over the covering ellipsoids using 
the polyhedral formulation easily finds a separating supset. 



relaxation, on the other hand, allows a lot of freedom in the optimization 
and separation can easily be established, see Figure A. 7 (right). 

We conclude this section by stressing that even if the 5-procedure is only 
sufficient when there is more than one quadratic constraint, adding new 
constraints can never make the inequalities harder to satisfy. On the con- 
trary, adding new terms may allow separations that would otherwise not be 
possible. 



A.6 Comments and References 



Computational Tests for Copositivity 

We have used the *S-procedure to enforce positivity of quadratic functions 
on polyhedra. It appears that this technique was first used by Shor in the 
late 80’s [187] and similar ideas appear in, for example, [162, 167]. 

Recently, Parrilo [ 1 54] has proposed a set of computational tests for copos- 
itivity based on the sums-of-squares decomposition of multivariable poly- 
nomials. These tests are formulated in terms of linear matrix inequalities 
and are strictly stronger than the simple decomposition results used here 
(for example, the approach can be used to establish copositivity of the ma- 
trix used in [48] to prove that the decomposition (A. 30) breaks down in R 5 ). 
These procedures could potentially be used to formulate stronger, but com- 
putationally more demanding, stability tests than the ones developed in this 
book. 
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B.1 Proofs from Chapter 2 



Proof B.l— P roposition 2.1 

Consider two polyhedral piecewise linear systems 



Si : 



x — XjX “T d i ~\~ BiU\ 
yi = CiX + Ci -p DiU\ 



S 2 : 



Z — X . j -p dj -p BjU 2 
V2 = CjZ + Cj + DjU2 



X G Xi i G I 
z G J G J 



where X* and denote polyhedral sets, represented as 



— {x | GiX -\- gi >: 0} Xj — | Gy 2 ; + S 0} 

When we interconnect the systems, we will obtain a partition with cells 
Xij = {(x, 2 ;) | x G Xi A z G Zj} 
which can be represented as 



X 



13 



x 

z 



Gi 0 
0 Gj 




i G /, j G J 



This gives rise to a partition with \I\x\J\ cells. Within each cell, the dynamics 
of each subsystem is linear and the dynamics of the interconnected system 
can be constructed using standard operations for linear systems. 

M. Johansson: Piecewise Linear Control Systems, LNCIS 284, pp. 171-184, 2003. 

© Springer- Verlag Berlin Heidelberg 2003 




Appendix B. Proofs 



For the series connection, we let u<i = y\ and obtain 



x 

z 



Ai 


0 


ai 




X 


BjC'i 




O/j + Bj Ci 




z 

_1 



Bi 

BjDi 



U\ 



X 



z 



e Xij 



V2 



D j Cj Cj* Cj | j Cj 



X 



z 



+ DjDiUi 



1 



x 



z 



e 



The parallel and feedback interconnections follow similarly. 



□ 



B.2 Proofs from Chapter 4 



Proof B. 2— Lemma 4.1 

The following formula for evaluating a piecewise smooth function W (t) will 
be useful. Let W(t) be a piecewise smooth function, and let t k denote the 
points of discontinuity of W(t). Then, 



W(t~) = W( 0) + 



±W(s)ds+ AW ( f k) 

k:tk<t 



(B.l) 



where A W(t k ) — W(t'l ) — W(tt_)- Obviously, if W(t) is non-increasing, then 
AW (t k ) <0 for all k and we have 

W(t~)<W( 0)+ / A W(s)ds . 

Jo ds 

Let us first consider the case a > 0. The inequalities (4.1) and (4.2) imply 
jV(t) <-^\\x{t)Y <-j^V{t) a.e. 

and multiplication with the positive function e 7t / ^ gives 

+ jvm = j t (^ t/P V{t)) < 0 a.e. 
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Letting W (t) = e 7 */ &V (t) and noting that W (t) is non-increasing we have 

ft j 

V(t~)e' 1t ^ <V( 0)+ / —V(s)e^ s ^ds. 

Jo ds 

The inequality (4.1) and the fact that V (t) is non-negative implies 

V(t~)e' 1t ^ < V(0) - 7 / V{s)e' 1s ^ds < V(0) 

P Jo 

Estimating V(t~) and F(0) using the lower and upper bounds in (4.2) gives 

i*wr < ^ix(o)ir e -^ 

a 

which establishes exponential convergence. 

Now, let the maximal a that satisfies (4.2) be negative. This implies that there 
is a time t = t 0 such that V(to) = a|x(t 0 )|| < 0 . Similarly to above, we have 

V(t~) < V{to) + f f-V(s) ds < V(to) + rr / V(a) da. 

Jto ds |a| J to 

Since V(t) is nonincreasing, we have V(t) < V(to) for t > to and 
V(t)< (l + (i-to)^|)v(io). 

Thus, V{t) — > — oo as t — > oo and by (4.2) it follows that ||x(t)|| p — > oc. 

□ 



Proof B. 3— Proposition 4.2 

Consider the Lyapunov function candidate V (x) = x T Px. Since P is posi- 
tive definite, V(x) can be bounded in the sense of Lemma 4.1, Equation (4.2). 
A solution to the strict inequalities (4.6) implies the existence of a 7 > 0 such 
that 

A? P -\~ P Ai -T 7 1 0 i — 1? • • • j !*• 

Now, consider the representation (4.5) of the differential inclusion (4.4). For 
the suggested Lyapunov function, V(x) = x T Px, we have 

7 L 

—V(x(t)) = '^2\ i {t)x{t) T {Aj P + PAi)x(t) < 

i= 1 

< XAOOlPOIl!) = -l\\ x( X)\l- 
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The desired result now follows from Lemma 4.1 by letting the vector norm 
| • | be the two-norm || • I 2 and by setting p = 2. □ 



Proof B. 4— Corollary 4.1 

Note that we can write the system (2.3) on the form (4.5) for almost all t by 
setting Xi(t) = 1 if x is in the interior of cell Xi and A i{t) = 0 otherwise. Since 
trajectories in the sense of Definition 2.1 do not remain on the boundary for 
any time interval, we do not need to define the dynamics on the cell bound- 
aries. Continuity of V(x) and x(t) implies that V(t) = V(x(t)) is continuous 
for all t. The result now follows from Lemma 4.1. 

□ 



Proof B. 5— Proposition 4.3 

If there exist Ri that solve the above inequalities, then for every P > 0 



0 < Tr [Pj2(RiAf + AiRi )] = J2 Tr l R i( A JP + PAi )] 

i i 



Hence, 0 < Tr [Ri(AjP+PAi)\ for some i, so AjP-\-PAi cannot be negative 
definite. □ 



Proof B. 6— Lemma 4.2 

Clearly, V{x) is piecewise quadratic. Continuity of V(x) follows from (4.14). 
The only obstacle in establishing (4.17) can occur around the origin. How- 
ever, the zero interpolation property guarantees that V (x) has no affine terms 
in regions that contain origin, i.e., 

V(x) = x T PiX = x T FjTFiX = x T FjTF iX = x T PiX x G Xi 

and the desired result follows. □ 



Proof B. 7— Lemma 4.3 

A solution to (4.23) guarantees that there exist > 0 such that 

z T (Pi - EjWiEi)z > tiZ T z 

for all z € R n+1 . In particular, for all x with x £ X it a solution implies that 

V(x) = x T PiX > x T EjWiEiX + CiX T x >0 x £ Xi 
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The last inequality follows from the fact that EiX y 0 for xgI* and that Wi 
has non-negative entries. Similarly a solution to (4.22) guarantees that 



V(x) = X T PiX > CiX T X 



x e Xi 



Thus, 



• f V(x) 
a = inf — jj 
xyo ||*|| 2 

is positive and finite, and satisfies V(x) > a||x|| 2 . 



□ 



Proof B. 8— Theorem 4.1 

Consider the Lyapunov function candidate V(t) = V (x(t)) defined by (4.24). 
Since trajectories x(t) in the sense of Definition 2.1 are continuous and piece- 
wise C 1 , Lemma 4.2 implies that V(t) constructed in this way is continuous 
and piecewise C 1 . Moreover, according to Lemma 4.2, there exists (3 > 0 such 
that the upper bound of (4.2) in Lemma 4.1 holds. A solution to the inequali- 
ties above implies that there exists a > 0 and 7 > 0 such that «!#(£) I 2 < V (t) 
and dV(t)/dt < — 'y ||o?(t) ||| - Hence, exponential convergence follows from 
Lemma 4.1 with || • || = || • I 2 and p = 2. □ 



Proof B. 9— Lemma 4.6 

We will first prove the following lemma. 

Lemma B.l 

The following statements are equivalent 

1. p T x > 0 for all x such that Ex y 0 and x ^ 0. 

2. p £ Int (JCe) where JCe = {y \ y = E T u , u y 0} 

3. There exists a vector u y 0 such thatp - E T u = 0 . 

□ 

Proof: Let 1C° E denote the polar cone of JCe, i.e., 

fC° E = {x | x T y < 0 , Vy G JC E } 
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and let B denote the ball B = {x | ||x|| = 1} Then, equivalence of Claim 1 and 
Claim 2 follows from 



1 p T x > 0 

rp X 

Op 1 — 

|| a; | 

O- p T x > 0 
<0 3e > 0, p T x > e 
■O p T x > e||x|| 
p T x > ey T x 
O- (p — ey) T x > 0 



Vx e -/C^\{0} 

V* € -/C^\{0} 

Vx € n B 

Vx g n b 

Vx G -JC° E 

Vx g -lC 0 E ,VyeB 

Vx g -^B,Vy g B 



Hence, by Farkas’ lemma, p — ey e ICe^V € B which implies 
o>p € Int(/C B ) o> 2. 

We prove equivalence of Claim 2 and Claim 3 in two steps. First 



2 ^ p - ty eJC E 



VyeB 



r i T 

Let 1=1 ... 1 . Then, the statement above implies 

pjT i 

= ^ p_e p 5 Tf € ^ 

pjT i 

O- 3m 0 > 0,P - e = E t u 0 

* p = e t ( U0 + e jWni) :=eTu 

where u >- 0. Hence Claim 2 implies Claim 3. Conversely, 



3 Op — E t u , w >- 0 
o3e > 0, u + ev >z 0 
<^36 !> 0,p -\- cE^ v G Ee 
Op + eE T v g 1C e 



\/veB 



and, if E has full column rank, 

Op G Int(/Cs) ^ 2 
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We can now proceed to prove Lemma 4.6. Note that Claim 2 implies 
Claim 1 trivially, since 



p T x = u t Ex > 0 
for u y 0 and all x with Ex ^ 0. 

Consider the converse statement. If E has full column rank, then 1=> 2 by 
Lemma B.l. When E does not have full column rank then we can, without 
loss of generality, assume that E is on the form 



E =g 




0 • 



where E + has full column rank. Now, Claim 1 implies that p must be on the 
form 



p = 



T 

pX 




Let x = |x + xoj • Then, we have 

> 0 Vx + with E+x+ > 0, 7^ 0 

and, by Lemma B.l, there exists >- 0 such that 

P+ - E\u + = 0 



Hence Claim 2 follows with u = 



X + U 0 



y 0, with arbitrary but uq >- 0. 



B.3 Proofs from Chapter 5 



Proof B. 10— Theorem 5.1 

Let i(t) be chosen so that x(t) e X if and let Pi = [I 0 ] T Pi [I 0] fori e /o. 
It then follows from the matrix inequalities in Theorem 5.1 that 



0 > 



PiAi + Af Pi + CfCi + EjUiEi 

B? Pi 



PiBi 

-7 2 I 



Vi el. 
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Multiplying from the left and right with [x T u] and removing the non- 
negative terms x T EjUiEix gives 

0 > 2 x T Pi(AiX + Biu) + x T CjCiX — 7 2 u T u = 

= + \\y\l -i 2 \\u\l 

Integration from 0 to oo gives the desired inequality. □ 



Proof B. 1 1 — Theorem 5.2 

It follows directly from the two inequalities that 

0 > P%Ai + AJ Pi + Qi + EffUiEi) i g I. 



Let i(t) be chosen so that x(t) £ Xi. Then, multiplying the above inequality 
from left and right by x and removing nonnegative terms gives 



0 > x T PiX ) + x(t) T Qix(t). 

Integration from t = 0 to t = oo gives the desired result. 



□ 



B.4 Proofs from Chapter 6 

Proof B. 12— Theorem 6.3 

It follows directly from the matrix inequalities in Theorem 6.3 that 

n ^ ~PiAi + AJ Pi EQi- EfUiEi PiBi] . 

Bf Pi Ri 

Multiplying from left and right by (x, u) and removing the nonnegative 
terms including Ui gives 

0 < 2 x T Pi(AiX + Biu) + x T QiX + u T RiU 

= (x T Pix) + x T QiX + u T RiU 

Integration from 0 to oo gives the desired result. □ 
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B.5 Proofs from Chapter 7 



Proof B. 13— Theorem 7.1 
Define the function 



V(x)=pfx xGlj, iel. 

Recall that by the construction of the matrices Fi for simplex partitions, the 
entries of t specify the values of V at the vertices. More precisely, we have 
V(0) = 0 and V{vk) = tk - The condtion 0 -< t -< 1 implies that 

ClIMloo < V(x) < C2 Halloo 

for some positive scalars ci, c 2 > 0. Denote the approximation error by 



<h(x) 



f(x) - AiX - <ii 

0 



x G Xi , i £ I. 



We have 

jV{x) = pj f(x) 

= pi(AiX + ai(x)) 

^ T - (5)1x100 -)- || ^ || oo * ||flj(x)|oo 

< — <5|x|oo < -8V(x)/c2. 

for some 5 > 0. Exponential stability follows. 



□ 



Proof B. 14— Theorem 7.2 

By a standard converse Lyapunov theorem, such as Theorem 3.12 in [116], 
there exists a C 1 Lyapunov function V (x) that satisfies 



ci\x\l < V(x) < c 2 \x\l 

f ( x ) < -C 3 \x[l 



dV_ 

dx 



< Ci\x\ 2 



(B.2) 

(B.3) 



for some positive constants c\ , <■■>, <’:>> and c\. The function V(x) can be ap- 
proximated by a function V on the form 

V(x) = x T PiX x £ Xi, i £ I 
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by letting Pi = FfTFi with Fi computed using Proposition A. 12 and 



Then, 



T = 



V{y x ) ... V{y p ) /2+ 



V(vi) ... V{y p ) 



T r 



1 ... 1 



/ 2 



= Tu := y(z/<) 



* = 



and y and <9y /<9x become arbitrarily accurate approximations of y and 
dV/dx as the partition is refined. Let 7 * be defined by the size of V. For 
sufficiently small approximation errors e* inequalities (B.2) and (B.3) imply 

ci\x\\/2 < x T PiX < c 2 \x\l 
x T (AJPi pPiAi)x < -2ei7<|x||2 



What remains is to find Ui and Wi with non-negative entries such that (7.6)- 
(7.9) hold. By the C 1 condition on / it can be assumed without restriction 
that V and y are quadratic and positive definite in a neighbourhood of the 
origin. Hence, Ui and Wi are not needed and can be put to zero. In the re- 
gions that do not contain the origin, V is linear, so Ui and Wi exist by Farkas’ 
lemma. □ 



B.6 Proofs from Chapter A 

Proof B. 15— Proposition A. 8 

We will first establish the following results. 

Lemma B.2 

Let x, 2 G M n and let x = Tz = Tz + 1 be an affine bijective map. If 

£(Si) = {z | z T §iZ < 1 } 

is the minimal volume ellipsoid that contains the set Z if then 
S{f T Sif) = {x | x T T T SiTx < 1} 
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is the minimum volume ellipsoid that contains the set 

Xi = {x = Tz \ z e Zi}. 



Proof: Since the mapping x = Tz, is affine and bijective, 

S(Si) DZi& £{T T SiT) D X t . 

Moreover, if 5(5*) describes an ellipsoid, it can be written as 
S(Si) = {z | (z - z 0 ) T P~ 1 (z - z 0 ) < 1} 
for some xq and some P = P T > 0. The volume of £(Si) is then 

vol(£(Si)) = Vdet PV n 

where V n denotes the volume of a unit sphere in R n . Let 



T = 



T 



r • 



□ 



Then, the volume of the transformed ellipsoid is 

vol(S(T T SiT)) = det T V det P • V n . 

Hence, for a given mapping T the volume is proportional to det P, and the 
circumscribing ellipsoid of minimal volume can be obtained by optimizing 
the volume of either 5(5*) or 5 (T t 5^T). 

Lemma B.3 

The minimum volume ellipsoid containing the standard simplex, 



n+ 1 

Xi = {x e R n+1 I X k > 0, = 1} 

k= 1 



is the ball 



n+l 

£* = {x \ x T x < 1, Xk = 1} (B.4) 

k= l 

□ 
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Proof See [98]. 



Proposition A.8 now follows trivially, since the map x z defined by 

z = Gix 

is a bijective affine map, that transforms an arbitrary simplex in R n into a 
regular simplex in the constraint hyperplane. Y2k=l z k = 1. By virtue of 
Lemma B.3, the circumscribing ellipsoid with minimal ellipsoid is then the 
ball (B.4). By virtue of Lemma B.2, the minimal ellipsoid that contains an 
arbitrary simplex in R n is given by 

£* = [x | x T GjGiX < l} . 

This concludes the proof. 

□ 



Proof B. 16— Proposition A. 10 

Since the minimum volume n-dimensional ellipsoid containing the unit cube 

Xi = {x e R n 1 1| 3 ? 1 00 < 1} 

is the ball {x e R n | x T x < n}, the results follows by direct application of 
Lemma B.2. □ 



Proof B. 1 7— Lemma A. 1 

Consider a cell Xi with cell identifier G- h partitioned as 



G% — 



Gi 

Hi 



0 

hi 



(B.5) 



First, consider i e Jo. By deleting all rows whose last entry is non-zero we 
find the cell identifer 



Ei = 



Gi 



0 



which clearly has the zero interpolation property. This proves the first claim. 
Now consider the second claim. The implication (A. 20) => (A.19) is imme- 
diate since Xi C {x \ EiX >z 0}. To prove the converse, (A.19) ==> (A.20), 
note that since 0 € Xi we must have hi >- 0 in (B.5). Hence, for every 
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x € {x | EiX y 0} there exists some e > 0 such that ex G Ij. Since V is 
quadratic we have 

V(x) > 0 ==> V (tx) >0 Vt 7^ 0 

In particular, since V (ex) > 0 we have V (x) > 0. This proves the second 
claim. For i e h we simply note that the added constraint [0 l]x > 0 is 
satisfied for all x e R n so Xi = {x | E % x y 0} and, hence, (A.19) (A.20). 

This concludes the proof. □ 



Proof B. 18— Proposition A. 14 

As discussed in the text, V (x) > 0 for xel*, i e I is equivalent to 

z[TiZi >0 % y 0 (B.6) 

where 2i = Gix are the Barycentric coordinates of x e Xi. Let O be a matrix 
where every entry is unity. Then, this inequality implies that there exists an 
ti > 0 such that 



z[ (Ti - 6i(I + 0 ))% >0 Zi y 0. (B.7) 

In fact, since z [(0 + I)zi > 0 on the domain of interest, e* can be taken as 
the minimal Ci that satisfies (B.7) on the compact set 

{ Zi | ||2i || 2 *= 1, % y 0} 

This implies that Ti — ^(7 + O) is copositive and, by Proposition A. 13, there 
exists Ui y 0 such that 

Ti — Ui — tiO — ( I > 0 



Let Wi = Ui + 6 iO. Then Wi >- 0, and 

Ti — Wi > til > 0 

Pre- and post-multiplication with the full rank matrix Ei and invoking the 
identity EjTiEi = FjTFi gives 

FjTFi-EjWiEi >0 

For i G Jo, we can follow the lines of Lemma A.l and establish a similar 
identity to (B.6) in terms of Zi = EiX (the partition coordinates for the non- 
zero vertices of Xi). The proof then follows similarly. 

□ 



183 




Appendix B. Proofs 



Proof B. 19— Proposition A. 15 

According to Proposition A. 8, the ellipsoidal bounding of a simplex that has 
minimal volume is given by 



Si = 



Onxn 9 

Ol xn 1 



- Ef Ei 



Let 0 = 



and thus 



'1 xn 



and 1 = 



L 1 xn 



. From the definition of Ei we have 



V E, = 0 



UiSi = mEi (IP - I)Ei = Ef mill 1 - I)E i 



which is on the form EjUiEi with = ui if i ^ j and 0 otherwise. This 
concludes the first part of the proof. For the second part of the proof, con- 
sider the simplex 



Xi = < x | x g co( 



for which we have 



Ek = 


"-1 

1 


-1 3 " 

0 -1 


, Si = 


"2 1-4" 

1 2 -4 




0 


1 -1 




-4 -4 10 



Let 



Pi = 



20 0 
0 1 
5 0 



5 

0 

-25 



Pre- and postmultip lying the LMI (A. 3 5) by z 



3 6 




we obtain 



-64 — 2 u if which is negative for all admissible values of ui (ui > 0). Hence, 
there is no solution to this LMI. However, for the formulation (A.36), it is 
straightforward to verify that the choice 



Ui = 



0 

20 

5 



20 

0 

20 



5 

20 

0 



solves the LMI. 



□ 
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